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Synopsis 

Name of student : Shiv Narayan Roll No. 9110470 

Degree for which submitted : Ph D Department: Electrical Engineering 

Thesis Title : PARAMETRIC SYNTHESIS OF AUTOPILOT 

FOR BANK-TO-TURN MISSILE 

Name of thesis supervisor : Dr K. E. Hole 

Month and year of thesis submission : July 1998. 

Two basic methods of controlling a missile to achieve the commanded acceleration 
are skid-to-tum (STT) and bank-to-tum (BTT) . In STT, the body acceleration is 
attained by permitting the missile to develop both an angle of attack and a sideslip 
angle . In contrast , a BTT missile should not have, ideally, any sideslip . To achieve 
the desired orientation, it is rolled (banked) such that the plane of maximum 
aerodynamic normal force is first oriented in the desired direction and then the force 
is controlled by adjusting the angle of attack. 

The function of an autopilot is to cause the missile to achieve the commanded 
acceleration as closely as possible , while maintaining closed loop stability and certain 
constraints on other operating variables. A BTT missile is asymmetric in configuration, 
having very strong roll-yaw coupling. Due to the asymmetric airframe, the missile 
has tendency for sideslip. If sideslip is allowed to develop, it would give rise, in 
conjunction with angle of attack , to large roll rates. Therefore, for a BTT missile, 
the -most significant constraint is on allowable sideslip. 
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Because of the large roll rates due to sideslip, the roll-yaw channels are 
strongly coupled and can be considered as a single system. The linearized pitch 
dynamics, on the other hand, can be treated separately from roll-yaw dynamics. 
Thus, on the whole, the autopilot design for a BTT missile consists, essentially, of 
the design of two autopilots : one for the pitch acceleration command tracking, known 
as longitudinal autopilot, and another for commanding the roll rate about the velocity 
vector of the missile, known as the lateral autopilot. 

There are basically two schemes of designing a linear time invariant (LTI) 
controller : the first one based on analytical methods and , the second one ba.sed 
on parameter optimization methods. The most common examples of analytical methods 
are linear quadratic Gaussian (LQG) based designs and optimal controller designs. 

These analytical methods minimize a simple objective function for finding a controller. 
The weight matrices are specified in the objective function. The major disadvantage 
of these methods is that the design specifications need to be translated into objective 
function and the weight matrices . In complex systems, this is a difficult problem . 

On the other hand, in the parameter optimization methods one starts with the 
controller structures . The next step is to form an objective function that represents 
design specifications. An objective function can be foitned by weighted sum . or 
maximum of various performance indices . Certain explicit constraints may also be 
added and the design is posed as nonlinear optimization problem which can be .solved 
using optimization methods. The main advantages of parameter optimization are : 

(i) A much greater range of objective functions and constraints can be allowed 
than those in analytical methods. 

(ii) If the designer is certain that some controller of the selected architecture 
will do the job, the parameters of the controller can be obtained using parameter 
optimization . 

(iii) If a controller has been designed using analytical methods, it can be 
improved using parameter optimization methods. 

(iv) The designer has the freedom as regards the structure of the controller. 
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However, the parameter optimization methods pose the problem that they do 
not guarantee global solution. Moreover, they require initial guess to start the 
optimization process. These drawbacks can be overcome by using recent technique 
of Simulated Annealing (SA), which provides global optimal solution. As regards the 
choice of the controller structure , the controller parametrizations based on the 
factorization approaches due to Youla et al and Desoer et al can be used. 

A fast maneuvering , high performance BTT missile has to function in an 
uncertain and changing environment. This demands certain level of time domain and 
frequency domain responses from the missile system. The parametric methods ot 
controller design provide a good framework for meeting the demanding requirements 
on the time domain and frequency domain performances. The time domain and 
frequency domain performance functions can be defined based on the design 
requirements and the design problem can be posed as an optimization problem that 
can be solved using optimization methods. Most of the autopilot designs for BTT 
missiles reported in the literature are based on the analytical methods. The autopilot 
design in this framework has not been reported . 

Therefore, the motivations behind the work reported in this thesis are : 

(i) To explore the possibility of parametric synthesis of longitudinal (pitch) 
autopilot and develop design algorithm using nonlinear optimization methods. 

(ii) To explore the possibility of parametric synthesis of longitudinal autopilot 
in two degrees of freedom (TDOF) configuration and develop the design algorithm 
using Simulated Annealing (SA) technique. 

(iii) To explore the possibility of parametric synthesis of lateral (roll-yaw) 
autopilot and develop the design algorithm using SA technique. 

(iv) To investigate the robustness of the designed autopilot. 

A brief description of the work reported in this thesis is given below . 
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The first chapter introduces the control problem for BTT missile, the autopilot 
design requirements, and presents a brief literature survey on the subject. The 
motivations behind the research work carried out in this thesis are given alongwith 
the brief description of the subsequent chapters. 

In the second chapter, the longitudinal autopilot has been designed based on 
the pole-zero assignment using nonlinear optimization methods. The control law is 
based on the Youla factorization in the ring of polynomials. The design problem has 
been posed as nonlinear optimization problem which has been solved using the direct 
search methods of Simplex search, Powell’s conjugate direction and Pattern search. 

In the third chapter, the longitudinal autopilot has been designed in TDOF 
configuration. The time domain and frequency domain performances have been achieved 
independently through the choice of two free parameters. These two free parameters 
are selected using closed loop model matching concept and frequency domain loop 
shaping /using constrained optimization with SA technique , for achieving the time 
domain and frequency domain performances, respectively. 

In the fourth chapter, the lateral (roll-yaw) autopilot has been designed for 
stability axis roll rate command tracking using constrained optimization with Simulated 
Annealing. The control law, which is based on Youla parametrization, facilitates 
input-output decoupling of roll and yaw channels. As a result of which body axis 
roll rate and yaw rate capabilities are achieved independently. 

In the fifth chapter, the robustness analysis has been carried out for neglected 
and mismodeled actuator dynamics and neglected sensor dynamics. 

In the sixth chapter, the thesis is concluded and the main findings and contributions 
are highlighted. 
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"skidding" motion to the missile; hence the name "skid-to-tum". In contrast, a 
BTT should not have, ideally any sideslip. To achieve the desired orientation, it 
is rolled (banked) so that the plane of maximum aerodynamic normal force is 
oriented in the desired direction and the force is controlled by adjusting the 
angle of attack. 

According to control of missile motion in pitch, yaw and roll plane, the autopilots 
are classified as longitudinal, lateral and roll autopilots respectively. 

The actual missile dynamics is nonlinear in nature and involves coupling among 
the pitch, yaw and roll axes. Moreover, the missile has to function in changing 
environments which gives rise to uncertainties in aerodynamic characteristics, mass 
and balance, wind, flexible body mode dynamics, actuator and sensor nonlinearities 
and noise. 

Because of the above characteristics of the missile plant, the possible approaches 
for the design of autopilot for the missile could be: 

(i) to design an autopilot for combined three axes nonlinear system to track 
acceleration commands in all the three channels (pitch, yaw and roll) and, simultaneously, 
ensuring stability of the overall system. 

(ii) to linearize the missile dynamics at an operating point and design linear 
time invariant (LTI) autopilot to track the commanded accelerations. 

When the design is based on linear time invariant (LTI) model, the dynamic 
interactions among pitch, yaw and roll axes could be ignored and the robust controller 
be designed to achieve the stability and performance requirements. In this approach 
coupling is taken as disturbance acting on the system and robust controller is designed 
for each channel to take care of the plant uncertainties and parameter variations. 

Another approach could be to take linearized plant as multi input and multi 
output (MIMO) system and design a linear time invariant controller for this plant 
for commands tracking, maintaining stability of the overall closed loop system. 

In an asymmetric bank-to-tum (BTT) missile configuration, as shown in Fig. 1.1 
[7], very strong roll-yaw coupling is present. In its operation the BTT missile is first 
banked so as to have plane of the maximum aerodynamic normal force in the direction 
of target and, then, develops an angle of attack to achieve desired acceleration. The 
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function of the autopilot is to cause the missile to achieve the commanded acceleration 
as closely as possible, while maintaining closed loop system stability and certain 
constraints on the other operating variables, and control signals. Due to asymmetric 
airframe, the missile has tendency for sideslip. If sideslip is allowed to develop, in 
conjunction with angle of attack it would give rise to large roll rates. Therefore, for 
a BTT missile, the most significant constraint is on the allowable sideslip. 

For the reason of large roll rates due to sideslip, the linearized pitch dynamics 
are separated from linear coupled roll-yaw dynamics. Therefore, the autopilot design 
for a kind of BTT missile shown in Fig 1.1, consists of the design of two autopilots: 
one for the pitch acceleration command tracking, known as longitudinal autopilot, 
and another for commanding the roll rate about the velocity vector of the missile, 
known as lateral autopilot. 


1.2 Literature Survey 

The autopilots for the missiles have been designed using classical, modem 
control, and adaptive control techniques. 

The vast majority of autopilots have been designed using classical approach 
[3,4]. In classical methods, the design is based on a linear missile model. The linear 
missile model is obtained by linearizing the dynamics of the missile at an operating 
point, and, the model thus obtained is assumed valid in the neighbourhood of the 
operating point. The design is performed by ignoring the dynamic interactions among 
the pitch, yaw and roll axes. This permits individual channel designs for each pitch, 
yaw and roll axes. 

Although, the classical designs are based on single input single output (SISO) 
control theory, their performance can be enhanced through multivariable analysis. 
Wise [4] used multivariable singular value techniques to maximize performance and 
stability robustness of an existing roll-yaw autopilot designed with classical method. 
Nesline et al [5] designed robust missile autopilot using combined optimal/classical 
approach. Modem control theory has been used in conjunction with classical control 
theory to achieve a robust design which deals with the uncertainties of the high 
frequency effects while giving optimal performance in the low frequency range. 
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William et al [6] designed an autopilot using modem state space approach. The 
design is based on a tenth order model of the missile dynamics which is decoupled 
into separate roll and pitch/yaw sub-systems by treating roll rate as an exogenous 
input to the pitch/yaw equations. The models are linearized over specified operating 
region and gains are scheduled on dynamic pressure and roll rate. Unmeasured states 
are estimated using a constant gain, reduced order Kalman filter, and the separation 
principle is invoked to compute controller gains. The desired roll channel time response 
is obtained using a pole placement controller design procedure. Robust pitch/yaw 
channel gains are determined using LQG theory with loop transfer recovery (LQG/LTR). 
Integral control action is obtained by treating the commanded accelerations from the 
guidance system as state variables. 

Other design methodologies using modem control theory are given in [7-15]. 
Wise [7] eliminated the steady state command errors, which are inherent in LQR 
error state equation formulations by incorporating integral control into LQG/LTR 
design. A trade off is established between stability robustness recovery (LTR) and 
sensor noise amplification, which limits the amount of recovery possible. 

Lin and Lee [10] used GSLQ (Generalized Singular Linear Quadratic) control 
to BTT autopilot design. Control constraints were imposed by pole placement techniques. 
Good stability margins were obtained at each of the output, and excellent command 
tracking was achieved in the presence of sinusoidal disturbances and roll, pitch and 
yaw couplings. 

Reichert [11] designed dynamically scheduling multiple linear time invariant 
(LTI) controllers, using //oo/p synthesis for systems with widely varying plant dynamics. 
Wise et al [12] and. Wise and Nguyen [13] designed autopilots using Hoo and 
projective controls, respectively, for normal acceleration command tracking. Projective 
controls is used to retain the eigen-structure of an H<x, optimal state feedback controller 
using a low order output feedback compensator, and retains the performance robustness 
and disturbance attenuation properties of the state feedback design. 

Sobel and Cloutier [14] applied eigen-structure assignment to the design of an 
autopilot for the extended medium range air to air technology (EMRAAT) missile. 
Eigen-structure assignment is a state variable multi input multi output (MIMO) design 
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method that allows placement of the system eigen-vectors as well as their corresponding 
eigen-values. In this design the eigen-structure assignment feedback gains are computed 
by choosing desired eigen-vectors based on mode decoupling and using the orthogonal 
projection solution to compute achievable eigen-vectors. 

Shamma and Cloutier [15] applied gain scheduling technique for a missile for 
longitudinal autopilot. This design methodology does not involve linearization about 
an operating point. Rather, the missile dynamics are brought to a quasilinear-parameter 
varying (LPV) form via a state transformation. Once in quasi-LPV form a robust 
controller using |i-synthesis is designed to achieve the angle of attack control via fm 
deflections. The final design is an inner/outer loop structure, with the angle of attack 
control being the outer loop. The effect of the inner loop is to linearize the missile 
dynamics in an approximate and band limited manner thereby leading to a simplified 
outer loop design with guaranteed inner loop robustness properties. 

Benshabat and Chait [16] applied quantitative feedback theory (QFT) to a class 
of missiles that exhibit unstable and non-minimum phase behaviour over a large flight 
envelope. This technique overcame the problems related with gain scheduling. The 
advantages in using QFT are that (1) it is possible to investigate the feasibility of 
fixed controllers for the whole flight envelope, and, (2) it reveals tradeoffs between 
the size of flight envelope, controller complexity and closed loop specifications. 

Lin et al [17] designed a robust autopilot for the HAVE DASH 11 missile 
system. The controllers are solved using the generalized Hamiltonian approach which 
unifies a class of robust control designs in the same framework. Among the autopilot 
designs are the linear quadratic Gaussian (LQG) autopilot, generalized singular linear 
quadratic (GSLQ) autopilot, and Hoo autopilot. All the designs are formulated in the 
canonical form and solved based on the generalized Hamiltonian formulation. Krause 
and Stein [18] proposed a general adaptive control structure within which an autopilot 
design was accomplished with tuned system performance and robustness guarantees. 
Kamen et al [19] applied an indirect adaptive control technique to a BTT missile 
autopilot design. 

Apkarian et al [69] applied linear parameter varying (LPV) techniques to the 
global control of the missile. This technique is the extension of the standard Hoo 
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synthesis technique. The proposed methodology produces an LPV controller that is 
automatically gain scheduled along the trajectories of the plant. Gratt et al [70] used 
feedback linearization for the autopilot design for the boost phase of the advance 
kinetic energy missile. Application of the feedback linearization methodology permitted 
an autopilot design without considering varying operating conditions. Wise et al [71] 
used nonlinear optimal control to design a pitch plane flight control system for 
a high angle of attack agile missile. The nonlinear Hoo control law was obtained by 
approximating the solution to the Hamilton-Jacobi-Isaacs equation. The solution to 
the partial differential equation was obtained by using the method of characteristics 
and was numerically approximated using successive approximations. 

Carter and Shamma [72] used linear parameter varying transformations to design 
a gain scheduled autopilot for a bank-to-tum missile. The gain scheduled design did 
not require linearizations about operating points. This framework was applied to the 
design of a coupled longitudinal/lateral bank-to-tum missile autopilot. The pitch and 
roll-yaw dynamics were separated and robust controllers were designed using ii-synthesis 
for both the pitch and roll-yaw channels. Fu et al [73] proposed an adaptive robust 
neural network based control approach for a bank-to-tum missile autopilot design. 

Zhu et al [74] designed a missile autopilot for the angle of attack and normal 
acceleration tracking using extended mean assignment (EMA) control technique for 
the linear time varying (LTV) systems. The EMA control technique is based on a 
new eigenvalue concept, called series-D (SD) eigenvalue for LTV systems. Geng et 
al [75] proposed a fuzzy neural network control system architecture called the fuzzy 
cerebellar model arithmetic computer (fuzzy CMAC) neural network. This fuzzy neural 
network architecture exploited a synergism between the original CMAC neural network 
and the theory of fuzzy logic controller. 


1.3 Motivation 

There are basically two schemes of design of a linear time invariant (LTI) 
controller which stabilizes a given LTI plant, and meet some design specifications. 
The first is based on analytical methods for determining a controller which is optimal 



INTRODUCTION 


9 


in some well defined sense, and the second is based on parameter optimization. The 
analytical methods have closed form solutions whereas the parameter optimization 
problems are usually solved by numerical methods. 

The first scheme is called the optimal controller design. The most common 
examples are linear quadratic Gaussian (LQG) based designs, and Hoa optimal 
controller designs. These analytical methods minimize a simple objective function for 
finding a controller. The weight matrices are, usually, specified in the objective 
function. The advantage of these optimal controllers is that they always find stabilizing 
controllers. The great disadvantage of these methods is that the design specifications 
are to be translated into objective function, and the weight matrices. In complex 
systems, it is difficult to fit all design specifications into, objective function and also 
the translation of the actual design specifications into specification of weight matrices 
is difficult [50,51]. Similar comments apply to the other optimal controllers also. For 
example, in Hoo optimal control method, the actual design specifications are to be 
translated into the weighting transfer matrices. 

On the other hand, the parameter optimization methods for LTI feedback design 
start with the controller structures. Controller structure means a system model whose 
parameter values can be adjusted. The next step in the design using parameter 
optimization methods is to select an objective function that represents design 
specifications. An objective function can be formed by weighted sum or maximum 
of various performance indices. The weights define the relative importance of the 
various performances of the system. Certain explicit constraints may also be added. 

After a controller structure, an objective function and some constraints have 
been specified, the design is formulated as nonlinear optimization problem. This 
problem can be solved using optimization methods. The following are the main 
advantages of the parameter optimization methods: 

(i) A wide range of objective functions and constraints can be allowed than 
those available in analytical methods. 

(ii) If the designer is certain that some controller of selected architecture will 
do the job, the parameters of the controller can be obtained using parameter optimization. 
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(iii) If a reasonable controller has been designed using analytical methods, it 
can be improved using parameter optimization methods. 

(iv) The designer has the freedom as regards the structure of the controller. 

However, the parameter optimization methods pose some problems : 

(a) Initial guess is required in optimization methods. For simple systems, the 
initial guess is not a big problem but for complex systems it is a drawback. 

(b) The nonlinear optimization problem may be nonconvex in nature. The 
general algorithms may provide only local optimal solution but not a globally optimal. 

These drawbacks can be overcome by using recent techniques of Simulated 
Annealing (SA) and Genetic Algorithms (GA) [38, 52, 63]. These techniques provide 
global optimal or near global optimal solution. 

The factorization approaches, initially presented by Youla et al [29, 32] and 
subsequently by Desoer et al [30, 44], provide simple parametrizations of controllers 
that stabilize a given plant. The controller parametrizations provide the choice of 
controller structure. The controller structure can be chosen through the choice of 
certain parametrization, and the controller design can easily be posed as optimization 
problem. 

The design specifications of a feedback control system can be expressed as the 
frequency and time domain specifications or in other words, design specifications 
require the shaping of several frequency and time domain responses. These responses 
can be expressed as hard bound constraints. So all the design specifications can be 
expressed as hard bound constraints on the frequency and time domain responses. 
The objective function can be formed as a weighted sum or the maximum of various 
frequency and time domain performance functions. 

From the literature survey, as in Section 1.2, it is clear that most of the autopilot 
designs for various kinds of missiles have used analytical methods. The BTT missile, 
shown in Fig. 1.1, is a fast maneuvering, high performance missile which is meant 
for air to ground applications. It has to function in uncertain and changing environments. 
This demands certain level of time domain and frequency domain responses from 
the missile system. The parametric methods of controller design with nonlinear 
optimization provide good framework for meeting -the demanding requirements on 



INTRODUCTION 


II 


time domain and frequency domain performances. Moreover, the roll and yaw channels 
in the above configuration of the missile are coupled and require to keep sideslip 
small to avoid large roll rates. These constraints can be easily incorporated in the 
frequency domain and time domain performance functions. Parametric synthesis 
methods for autopilot design in this framework have not been reported in the literature. 

Therefore, the motivations behind the work carried out in this thesis are ; 

(1) to explore the possibility of parametric synthesis of longitudinal (pitch) 
autopilot and develop design algorithm using nonlinear optimization methods. 

(2) to explore the possibility of parametric synthesis of longitudinal autopilot 
and develop design algorithm in two degrees of freedom (TDOF) configuration using 
Simulated Annealing (SA) technique. 

(3) to explore the possibility of parametric synthesis of lateral (roll-yaw) 
autopilot and develop the design algorithm using SA technique. 

(4) to investigate the robustness of the designed autopilot. 


1.4 Thesis Organisation 

The first chapter introduces the control problem and autopilot design requirements 
for a BTT missile, presents brief literature survey on the subject, sets the motivation 
behind the research work carried out in this thesis. 

In second chapter, a longitudinal autopilot has been designed based on pole-zero 
assignment using nonlinear optimization methods. The control law is based on the 
Youla parametrization in the ring of polynomials. The design problem has been posed 
as nonlinear optimization problem which has been solved using direct search methods 
of Simplex search, Powell’s conjugate direction and Pattern search . The performance 
function includes the terms corresponding to the command tracking, overshoot, 
undershoot, gain and phase margins. 

In third chapter, the longitudinal autopilot has been designed in two degrees 
of freedom (TDOF) configuration. The time and frequency domain performances have 
been achieved independently through the choice of two free parameters. These two 
free parameters are selected using closed loop model matching concept and the 
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frequency domain loop shaping using constrained optimization with Simulated Annealing 
technique, for achieving time domain and frequency domain responses, respectively. 

In fourth chapter, the lateral (roll-yaw) autopilot has been designed for the 
stability axis roll rate command tracking using constrained optimization with simulated 
Annealing. The control law which is based on Youla parametrization facilitates 
input-output decoupling of roll and yaw channels as a result of which the body axis 
roll rate and yaw rate capabilities are achieved independently. 

In fifth chapter, the robustness analysis has been carried out for neglected and 
mismodeled actuator dynamics and neglected sensor dynamics using singular values. 

In sixth chapter, the thesis is concluded and the main findings and contributions 
are highlighted. 



Chapter 2 


Pitch Autopilot Synthesis Based on 
Pole-Zero Assignment with Nonlinear 
Optimization 

2.1 Introduction 

In control systems, the pole dominates the transient response and system stability 
[20-30], and, the zero of a system plays important role in the interaction between 
the system and environment [24-26]. So there have been many studies on the pole 
and zero assignment methods. Another strategy for control system design is the 
pole-zero assignment to minimize the control input signal of the system [25]. In 
optimal control theory [28] , specially LQR and LQG, control input signal is minimized 
in the performance index. 

One of the important results in the design of control systems is the parametrization 
of all stabilizing controllers in the unity feedback configuration [29,30]. 

With the development of modem computing facilities, these algebraic methods 
are advantageous in the following sense : 

(i) The designer’s level of subjective judgement is minimised. 

(ii) Design goals, performance as well as stability, can straight way be incorporated 
in the performance index and the problem can be cast as parameter optimization 
problem which can be solved using optimization techniques. Thus parametrized 
controller can be constructed in a systematic way to assure stability as well as other 
design goals. 
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(iii) The algebraic methods permit assignment of poles and zeros through the 
selection of polynomials. 

In this chapter, a parametrized longitudinal autopilot is designed using pole-zero 
assignment. The pitch dynamics form a single input multi output (SIMO) model. In 
the missile configuration, Fig.2.1, the rate feedback forms the inner loop whereas 
the acceleration feedback forms the outer loop. In this methodology one loop is 
designed at a time. While the inner pitch rate loop provides stabilization, damping 
and minimum input control, the outer acceleration loop is designed for acceleration 
command tracking with specified gain and phase margins. 

The pole-zero assignment control law is developed based on the method by 
Tu [27], which uses Bongiomo and Youla’s factorization approach [29] in the ring 
of polynomials. 

The overall closed loop transfer function places a pair of dominant complex 
conjugate poles to achieve performance measures e.g. rise time, percentage overshoot 
and settling time etc. The remaining system poles and zeros are then selected via 
nonlinear optimization to give desired gain and phase margins. The performance index 
includes the terms corresponding to command tracking error, overshoot, undershoot, 
gain and phase margins. As this performance index is nondifferentiable, the gradient 
based optimization methods would provide poor results. Therefore, direct search 
methods of Simplex search, Powell’s conjugate direction and Pattern search have been 
employed for optimization and their performance has been compared. 

The structure of the chapter is as follows ; in Section 2.2, the missile pitch 
dynamics is described. In section 2.3, the pole-zero assignment control law is discussed. 
In Section 2.4, the minimization of control input is described. In Section 2.5, the 
design algorithm is given and in Section 2.6, the design algorithm is implemented 
for the design of longitudinal autopilot. In Section 2.7, the design results are given 
and are compared with the results of other design methodologies. In Section 2.8, 
the chapter is concluded. 
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Figure 2.1 : Longitudinal Autopilot 

2.2 Missile Pitch Dynamics 

The details of the missile pitch dynamics, taken from [7,13], are described in 
Appendix 1. The equations governing missile pitch dynamics are rewritten as 
follows : 

d = Za-a + Z5-8e + 9 
q = M^-a + Afg-S^ 

5 = - 2 ^ CO 5^ - co^ (5 - 5^ ) 

where a,q ,he and (O are angle of attack, pitch rate, elevator fin deflection, 
fin actuator damping and natural frequency, respectively. The measured variables are 

normal acceleration, Az = VZa • cl + VZSe • 5e (ft/s ) and the pitch rate q. The scalar 
control input, u = hec{rad) is the elevator fin angle command. V is the missile 
velocity. 

The transfer function matrix from the elevator fin command, 5ec, to normal 

acceleration, Az, and pitch rate, q, is given by : 

' Az' (sftv (Z^-s'^ + ZaMh-Z^Ma) 

Sec {s^-Zas-Ma)(s^ + 2q(iis + d) 

Gp(s) = = . 

oft (Md- s+MaZs-MsZa) 

(ft-Zas-Ma)(s^ + 2cO)s + (£)^) 
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Figure 2.2 : Unity Feedback Configuration 


2.3 Pole-Zero Assignment Control Law 

Consider the unity feedback configuration shown in Fig.2.2 
Let the plant be 

where A(s) and B(s) are coprime polynomials. 

The reference input is 


m = 




M{s) ' 

where M(s) and N(s) are coprime polynomials. 
The sensitivity function is given as 
Sis) = 


( 2 . 1 ) 


( 2 . 2 ) 


(2.3) 


1 + Pis) Cis) 
and, the error is given by 

Eis) = Sis) Ris) (2.4) 

Let A(s), B(s) and M(s) be factorized as 

/4(j) = A+(i) • A-is) (2.5) 

Bis) = B+is) • B-is) (2.6) 

Mis) = M+is) ■ M-is) (2.7) 

where A+ is), B+ is) and M+ is) have all their zeros in Reis) > 0 while 
A- (s), B- is) and M- {s) have all their zeros in Reis) < 0. 

For a reference signal tracking and desired pole assignment, the sensitivity 
function must be of the form [27] : 

Wis) M+is) 


Sis) = 


Gis) 


( 2 . 8 ) 
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where G{s) is Hurwitz polynomial with desired closed loop poles and W( 5 ) is 
undetermined polynomial which should satisfy the internal stability constraints. 

2.3.1 Definition [27] 

The sensitivity function S(s) is said to be internally stable if the closed loop 
system of Fig. 2. 2 is asymptotically stable for some choice of C(s) ie. there is no 
pole-zero concellation between C(s) and P(s) in Re[s] > 0. 

The following Lemma gives conditions which guarantee internal stability of the 
sensitivity function. 


2.3.2 Lemma [27] 

The sensitivity function S(s) is internally stable if and only if all the following 
conditions hold: 

(i) S(s) is analytic in Re[s] > 0 

(ii) The numerator polynomial of S(s) is divisible by A+(5). 

(iii) The numerator polynomial of l-S(s) is divisible by B+{s). 

Equation (2.8) and condition (ii) of the Lemma imply that the numerator of 

S(s) must contain the least common multiple (LCM) of A+( 5 ) and M+(s). ie. 

W(^)M+(^) L{s)Z{s) (2.9) 

GW = GW ■ 

where Z(s) is the least common multiple of A+{s) and M+(.s) while L(s) is 
undetermined polynomial and is used for loop shaping. 

To satisfy the requirements of causality, S(s) must be proper i.e. 

deg ids)) > deg ms)) + deg (Z(^)) (2. 10) 


From the equation (2.9), we have 

.f,.G{s)-L(s)Z(s) 

T{s)-l-S{s)- 


( 2 . 11 ) 


From condition (iii) of the Lemma, we have 
His) = Gis)-L{s)Zis) 


= B^is)Fis) 


( 2 . 12 ) 
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where, F(s) is undetermined polynomial, the choice of which facilitates 
incorporation of additional zeros. 


2.3.4 Theorem [27] 

The solution of L(s) in equation (2.12) exists if and only if M+{s) is coprime 
with B+ {$). 


The controller, thus, can be derived from the equation (2.3) as 

c(.)= 

P{s)S{s) 

_ A{s)BJ,s)F{s) 

~ B{s)L{s)Zis) 

The output is given by 
y(5) = (l-5(i))-i?(j) 

_ BJ,s)F{s)N{s) 

G{s) M{s) 

The closed loop transfer function is given by 
B^{s) F{s) 


ns) = 


Gis) 


The loop transfer function, Li(s), at the input of plant, is given by 
Li(s) = P(s) C(s) = 


(2.13) 


(2.14) 


(2.15) 


(2.16) 


L(s) Z(j) 

The desired loop shaping and command tracking can be achieved by proper 
choice of F(s), Li(s) and Gis). 


2.4 Minimisation of Control Input ' 

The control input signal, u(t), is minimized to save the control energy and to 
avoid saturation problems. 

The L2 norm of u(t) is defined as [32] 

^ ^I/2 

1|m 1|2 = j u\t) dt 
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If u(t) is bounded for all values of r > 0 ; is zero for r < 0 and approaches 

zero at least as fast as e ^ ^ as t approaches infinity, where 6 is a small positive 

constant, then applying the Parseval theorem [39], L 2 norm of u(t) becomes 

< > 1/2 

r 

Il“ll 2 = j u{t) dt 
0 

\ J 

( > 1/2 

= I U{s)U{rs)cb 

2nj j 

-y 00 

> y 

(2.17) 

where U(s) is the Laplace transform of u(t), i.e U(s) is stable with poles in 
Re[s] <0. 


2.5 Design Algorithm 

The design algorithm is aimed at achieving closed loop pole-zero configurations 
for inner pitch rate loop and outer acceleration loop for the specified performance. 
The design is achieved in two parts: in Algorithm I, the inner pitch rate loop and 
in Algorithm 11, the acceleration command tracking loop is designed. 

In the design methodology given above, once the overall transfer function is 
chosen, the rest of the design for parametric controller is straight forward. The crux 
of the design is how to choose transfer functions for inner pitch rate loop and outer 
acceleration loop. The choices can be based on the minimization of the integral of 
time multiplied by absolute error (ITAE), quadratic performance index, model matching 
or computer simulations. 

Some constraints are imposed on the bandwidth of the overall system or on 
the actuating signal. For the inner pitch rate loop, since our concern is with stabilization, 
damping and control input minimization, the proper choice would be to choose 
transfer function based on ITAE and impose constraint on u(t). 
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Overall closed loop transfer function for the acceleration command tracking 
loop is obtained by model matching through nonlinear optimization and constraints 
are imposed on maximum overshoot, maximum undershoot, gain and phase margins. 

Let the plant transfer functions for the inner and outer loops be denoted by 
Pl(s) and P 2 {s), respectively. Ti(s) and T2{s) denote closed loop transfer functions 
for inner and outer loops, respectively. 

The L2 norm of the control input, u, to the actuator due to unit step input is 
minimized. This saves the control energy and avoids saturation problems [27]. 

Now design algorithms for the pitch rate loop and acceleration command tracking 
loop are given as Algorithm I and Algorithm 11, respectively. 


2.5.1 Algorithm I : (for pitch rate loop) 

Step 1 

Optimal ITAE, expressed as polynomial in coo is selected from the standard 
forms [33] for U(s) to be strictly proper. The standard form of ITAE, however, is 
not available for all plant transfer functions. In such cases the idea of [34] may be 
used and a digital computer simulation may be used to obtain its optimal transfer 
function under the constraint on u(t). 

Step 2 

By computer simulation, coo is obtained to give |u(t)|<1.0 for all t > 0. 

Step 3 

L2 norm of u(t) i.e. || u(t) 1|2 is minimized using nonlinear programming under 
the constraint of equation (2.12) thus assigning the zeros of the system suitably. 

Step 4 

After ascertaining closed loop transfer function Ti{s) and hence L(s), F(s) and 
G(s) from equations (2.11), and (2.12), the controller Ci(s) is parametrized using 
equation (2.13). 
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2.5.2 Algorithm II : (for acceleration command tracking loop) 
Step 1 

Choose the second order reference linear model whose time response behaviour 
is consistent with the specified performance measures of the missile system. 

Step 2 

Degree of the numerator and denominator polynomials, n(s) and m(s), of the 
closed loop transfer function T 2 (s) is selected such that : 

{Np — Wj)cIosed loop system ^ (A^p ~ N^optn loop system (2. 1 8) 

where Np and denote number of poles and zeros, respectively, i.e. 

nis) = (5 + ai) (j + <X2) ..... {s + aNz) 
mis) = (5 + pi) is + P 2 ) .... is + PATp) 

An upper limit is placed on the selection of a/ and P; to avoid numerical 
overflow during the optimization process and also because highly damped poles tend 
to affect the system dynamics adversely [35]. 

Step 3 

Search for at and P/. 

The aj and P/ are design variables denoted by the vector X. They are selected 
via nonlinear optimization techniques. 

The design variables are selected to minimize the following performance index 

^ “t 2 r ^ 

~ 2) \_y^ J ^ f Jmax ~ Jraax 1 + < 3^niin 

+ (ZLi (/'cOc) +180“- Qspecf + (20 logio I Li (/-CDg)! - AO^ ] (2. 19) 

where X is the penalty factor; y it) is the desired time response ; <ymin> is 
zero if ymin^O or otherwise it equals (ymin) • 

The first term in the penalty term provides a measure of the error for the 
percentage overshoot: 
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The third term in the penalty term accounts for phase margin (>Qspec)- The 
phase margin specification is satisfied when the phase angle of Lj(/C0c) is greater 
than or equal to (180° - G^pec) at frequency cOc when the magnitude of L,(/C0c) is 
equal to 0 dB. Liijo) is the loop tranfer function with loop broken at the input of 
actuator of the missile. 

The last term in the penality term imposes the constraint on the gain margin 
(> NdB) 

Step 4 

After ascertaining closed loop transfer function, Tiis) and hence L(s}, F(s) and 
G(s) from the equations (2.11) and (2.12), the parametric controller C2(5) is obtained 
using equation (2.13). 


2.6 Design Example 

Design example has the model described in Section 2.2 of this chapter. The 
flight conditions studied here are the same as considered in [13], represented by a 

= 16®, Mach 0.8 and an altitude of 4000 ft. The nominal values of the dimensional 

aerodynamic stability derivatives are : 

Za = -1.3046 (1/s) ; Z5 = - 0.2142 (1/s) 

Ma = 47.7109 (1/s^); and M 5 = - 104.8346 (1/s^) 

The remaining system parameters are ; 

V = 886.78 (ft/s); ^ = 0.6, (0 = 113.0 (rad/s) 

The transfer function matrix, Gp (s), is given by 

■ - 2.4254e6 (5 + 26.196) (^- 26.196) 

Sec {s + 67.8 + 90.4i) {s + 67.8 - 90.4/) {s + 7.5903) (s - 6.2857) 

Gp(s) = 

-1.3386e6 (5 +1.4021) 

Sec {s + 67.8 + 90.4/) (5 + 67.8 - 90.4/) {s + 7.5903) {s - 6.2857) 


( 2 . 20 ) 

The transfer function from fin command 5ec to Az has a RHP zero and the plant 
has also a RHP pole. Thus the missile pitch plant is unstable and non-minimum 
phase. 
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Typical design specifications for such a missile can be stated as : 

Design a parametric longitudinal autopilot to track the commanded acceleration 
maneuvers with a rise time (63%) less than 0.3 sec, overshoot less than 5%, initial 
undershoot less than 2.5%; The controller must provide gain and phase margins of 
(-6, +10 )dB and ± 60 degrees, respectively. The controller bandwidth is to be restricted 
by high frequency unmodeled dynamics constraints (eg. ummodeled flexible modes 
and actuator nonlinearities). 

The controller is synthesized using Algorithms I and II of the design algorithm 
discussed in Section 2.5. 

2.6.1 Controller Synthesis for Pitch Rate Loop 

Algorithm I is used for designing pitch rate loop. 

For the inner pitch rate loop, the plant transfer function is 

- 1.3386e6(j+ 1.4021) 

~ {s + 67.8 + 90.4i) (j + 67.8 - 90.4!) (s + 7.5903) (s - 6.2857) 

The controller Ci(s') is designed with || zi(?) ||2< 1-0, for all t > 0. 

From Fig. 2.2 the control input U(s), can be written in terms of sensitivity 
function S(s) as 

U{s) = C{s) S{s) R{s). 

Substituting the values of R(s), S(s) and C(s) from equations (2.2), (2.9) and (2.13), 
it reduces to 

^ A{s) F{s) N{s) (2.21) 

B-(s) Gis) M{s) 

Factorizing the plant transfer function Pi{s) according to equations (2.5), (2.6) and 

(2.7), we get 

B+is) = 1.0 

B-(s) = - 1.3386e6(s+1.4021) 

A+(s) = (s-6.2857) 

A-(s) = (s+67.8+90.4i) (s+67.8- 90.4i) (s+7.5903) 

Z(s) = s(s -6.2857) 

M+{s) = ^ 
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M-(s) = 1. 


The U(s) given by equation (2.21), fully contains the closed loop transfer 


function given by equation (2.15). Substituting the above factors into equation (2.21), 
it is clear that, deg(G(s))-deg(F(s)) should be at least 4 for the U(s) to be strictly 


proper. Therefore, with one zero, the optimal ITAE is taken of the form 

/: (i + Y) 


T\{s) = 


( 2 . 22 ) 


{s^ + 2.8 coo / + 5co§ + 5.5 (o3 + 3.4 co^ ^ + coB) 


The standard form for fifth order ITAE as given above is taken from [33]. To 
satisfy the equation (2.12), the L(s) and F(s) are chosen of the following forms: 
L{s) = 5^ + hP" + I2S + I3 
F{s) = fis+fi. 

Then H(s), given by equation (2.12) is 
His) = Gis) - u,s) Zis) 

= (/ + 2.8 cOq / + 5 cDq / + 5.5 coj s'^ + 3.4 (oj 5 + o) J - (/ - 6.28575) ( 5 ^ + + l^s + 1^) 

= B_^is) Fis) 

= Fis) 

= f,s^f2 
This gives. 


2.8 coo = - 6.2857 

5 a)§ = /2 - 6.2857 h 
5.5 (oB = Z 3 - 6.2857 h 
3.4o^t + 6.2S51 13 = fi 
(3^^ = fl 


(2.23) 


The control input to the plant turns out to be 

(s + 67.8 + 90.40 is + 67.8 - 90.40 is + 7.5903) (5 - 6.2857) if s +f ) 

-1.3386e6 5 (5 + 1 .4021) (5^ + 2.8cOqj'^ + 5 coJ 5^ + 5.5 s^ + 3.4 ©Jj + ©J ) ( 2 . 24 ) 

/ 

Now the problem is formulated as nonlinear optimization problem to minimize 
the L 2 -norm of u(t) i.e. 11(/(5)||2 under the constraints of equations (2.23). To obtain 
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the inital guess for coo, Y= 1.4021 is chosen which is existing zero of the plant. 
Computer simulation is used to find out the initial guess on coo to get 

1 I «( 0||2 < 1 . 0 . 

Simplex search with multiple restarts with X = 50.0, has the solution vector 
(coo , /l , /2 , /3 ,/l ,f2 } as {9.9988, 34.42, 716.245, 10000.2, 96843.0, 9.994e4}. This 
gives : 


Txis) = 


96843(5+1.032) 

(5^ + 27.997 + 499.88 + 5498 + 33984 i + 99940) 


and 


Cl ( 5 ) = 


-7.2346 e- 2 (5 + 7.5903) {s + 67.8 + 90.4/) {s + 67.8 - 9Q.Ai) 

(5+ 1.032) 

s{s + 1.4021) + 34.43 + 716.2455 + 10000.2) 


(2.25) 


(2.26) 


2.6.2. Controller Synthesis for Outer Acceleration Loop 


For the outer acceleration loop, the plant transfer function becomes 


P 2 {s) = -^-TUs) 
A-z/^ec 


Tl(s) 


q/^ec 

V J 

1.7547e5 (5 + 26.196) (5 - 26.196) (5 + 1.032) 

(5 + 1.4021) (5^ + 27 .9975"^ + 499.885^ + 54985^ + 339845 + 99940) 


(2.27) 

(2.28) 
(2.29) 


As the RHP zero is beyond the bandwidth of the system, there is no need to 
assign the zeros. However, LHP zeros can be assigned to new locations to get 
desired overshoot and settling time. 

The second order linear system which meets the time response specifications 
(desired), is given by 

(/+ 105 + 50) 

The zeros of the closed loop system are maintained at ai , a 2 = ± 26.196. 
Additional minimum four number of poles are required to satisfy the criterion. 
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(Np — yj)closed loop system ^ {Np - Nz)open loop system (2.30) 

These additional four poles are placed at —pi, -p 2 , —ps, and -pA to give : 

m{s) = + 10s + 50) {s+pi){s+ p 2 ) (s + ps) (s+pa) (2.31) 

The optimum set of intdependent variables X = \pi ,p 2 ,p 3 ,p 4 | is searched such 
that the performance index J given by equation (2.19) is minimum. 

' As J is non differentiable the gradient based optimization methods will give poor 
results. The direct search methods for nonlinear optimization are the ones which are 
suitable. For optimization of the performance index J, three direct search methods 
namely, Simplex search, Powell’s congugate direction and Pattern search have been 
applied and their performances have been compared. 

The Simplex method is due to Nelder and Mead [35]. This method is very 
robust but not very efficient in terms of the number of function evaluations needed 
for function minimization. The Powell’s conjugate direction method [37] is faster 
than simplex search. In this method the basic idea is to create a set of N linearly 
independent search directions, starting each time from the previous best point. 
Hooke-Jeeves Pattern search method [38] works by creating a set of search directions, 
iteratively. 

Solutions obtained using these techniques with different values of X and starting 
solutions are shown in the Table 2.1. The results are fairly insensitive to the type 
of optimization technique employed as well as to the starting solution. 

The optimal solution vector {pi ,p2 ,p3,p4} using Simplex search with multiple 
restarts (X = 10) is obtained as {41.223, 41.243, 41.236, 41.214}. 


This results in the closed loop transfer function 


-2.1079e5ri+1.4465e8 

Tiis) = 

(rf+ 10^ + SO) (s + 41.223) {s + 41.243) (s + 41.236) (s + 41.214) 


(2.32) 


The corresponding controller is 

^ -1.2013 {s + 1.4021) {s - 26.196) + 27.99?/ + 499.88/ + 5498/ + 33984^ + 99940 

■j(y + 1 .032) (/ l.T508e2 / + 1.192e4 / + 3.9161e5 / + 6.43419e6 5 + 4.3067e7) (2 33) 
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Table 2.1 




Variation of solution 

vector X = 

» ^2’ Py respect to starting solution and penality 

factor (X). 

X 

Starting 

Simplex Search 

PowelFs conjugate Direction 

Pattern Search 

Solution 

Variable 

J 

Variable 

J 

Variable 

J 

1.0 

30.0 

41.532 


45.778 


41.326 



30.0 

40.885 

2.471 

45.893 

2.192 

41.168 

2.496 


30,0 

41.087 

45.491 

38.116 


30.0 

41.165 


41.727 


40.616 



50.0 

41.500 


45.148 


41.502 


1.0 

50.0 

41.539 

2.464 

44.619 

2.190 

41.681 

2.465 

50.0 

41.451 

44.794 

40.980 


50.0 

41.497 


44.963 


41.550 



10.0 

41.528 


45.323 


41.632 


1.0 

30.0 

41.516 

2.4635 

44.631 

2.190 

41.568 

2.479 

50.0 

41.495 

44.337 

40.630 


100.0 

41.511 


45.135 


39.678 


5.0 

30.0 

41.339 


42.330 


41.832 


30.0 

41,164 

10.687 

42.637 

10.428 

41.681 

10.731 


30.0 

41.298 

42.032 

41.916 


30.0 

41.301 


41.152 


41.682 



50.0 

41.261 


42.010 


41.326 


5.0 

50.0 

41.329 

10.653 

41.724 

10.426 

41.613 

10.826 

50.0 

41.265 

41.723 

43.816 


50.0 

41.321 


42.317 


43.126 



10,0 

41.336 


42.213 


43.629 


5.0 

30.0 

41.303 

10.686 

42.192 

10.426 

41.668 

10.7 

50.0 

41.296 

41.965 

41.239 


100.0 

41.276 


42.018 


41.199 


10.0 

30.0 

41.223 


41.683 


41.261 


30.0 

30.0 

41.243 

41.236 

20.96 

41.777 

41.500 

20.703 

41.363 

41.268 

20.96 


30.0 

41.214 


40.934 


40.932 



50.0 

41.286 


42.105 


41-218 


10.0 

50.0 

41.304 

20.96 

41.585 

20.702 

41.206 

20.958 

50.0 

41.200 

41.442 

41.266 


50.0 

41.242 


41.743 


41.267 



10.0 

41.267 


41.807 


41.268 


10.0 

30.0 

50.0 

41.269 

41.270 

20.95 

41.633 

41.607 

20.701 

41.206 

41.266 

20.958 


100.0 

41.267 


41.344 


41.267 
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Simplex Search 


Variable 


Powell's conjugate direction 
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3 


Pattern Search 
Variable J 


Starting 

Solution 


30.0 

30.0 

30.0 

30.0 

50.0 
50.0 
50.0 

50.0 

10.0 

30.0 

50.0 
100.0 

30.0 
30.0 
30.0 

30.0 

50.0 
50.0 
50.0 

50.0 

10.0 

30.0 

50.0 
100.0 

30.0 
30.0 
30.0 

30.0 

50.0 
50.0 
50.0 

50.0 

10.0 

30.0 

50.0 
100.0 


41.254 

41.256 
41.263 
41.239 

41.257 

41.254 
41.252 

41.255 

41.266 
41.268 
41.270 

41.267 

41.237 
41.139 
41.176 
41.298 

41.242 
41.226 
41.181 
41.255 

41.249 

41.249 

41.246 

41.247 

41.236 

41.238 
41.208 
41.766 

41.249 

41.243 

41.268 
41.273 

41.216 

41.328 

41.242 

41.532 


41.507 


41.506 


41.506 


103.144 


103.143 


103.144 


205.879 


205.871 


205.873 


41.602 
41.691 
41.382 
41.009 

41.470 

41.407 

41.306 

41.432 

41.375 

41.482 

41.603 
41.414 

41.352 

41.490 

41.287 

40.868 

41.412 

41.203 

41.256 

41.179 

41.273 

41.368 

41.523 

41.251 

41.295 

41.505 

41.191 

40.755 

41.222 

41.056 

40.997 

41.253 

41.320 

41.252 
41.351 
41.307 


41.249 


41.248 


41.247 


102.888 


102.885 


102.885 


205.624 


205.614 


205.612 


41.256 

40.682 

40.732 

41.558 

41.256 

41.367 

40.676 

40.326 

41.326 
41.236 

41.266 
41.260 

41.268 

40.263 
41.319 
41.238 

41.249 

41.321 

41.386 

41.256 

41.268 

41.232 

41.169 

40.681 

41.257 

41.263 
41.256 
40.931 

41.251 

41.267 
41.259 
40.328 

40.691 

40.326 

41.369 

41.245 


41.507 


41.523 


41.506 


103.146 


103.145 


103.145 


205.908 


205.908 


205.924 
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Figure Response 

2.7 Results and Discussion 


Achieved design specifications with above solution vector are: 


Gain Margin 

= (-10.88, ~) dB 

Phase Margin 

= ± 63.95 Deg 

Loop gain cross over frequency 

= 6.9 Hz. 

Rise time (63%) 

= 0.32 sec 

Maximum overshoot 

= 4.06% 

Initial undershoot 

= 1.53% 

Settling time (95%) 

= 0.52 sec. 

Maximum fin deflection 

= 0.945 Deg/G 

Maximum fin rate 

= 5.96 Deg/s/G 


Unit step responses of the model and the present design methodology are shown 
in Fig.2.3. Nyquist plot of loop transfer function is shown in Fig.2.4. The elevator 
fin deflection is shown in Fig.2.5. 

The time and frequency response specifications achieved in the present design 
are compared with the other autopilot designs listed in Table 2.2, which is drawn 
from [13], having the same missile model. Design specifications achieved in the 
present design are comparable with any other design methodology. 
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Figure 2.4 ; Nycjuist Plot 
Table 2.2 


Autopilot Design 


* 

63% 

Rise 

Time 

95% 

Settling 

Time 

% 

Initial 

Under 

shoot 

% 

Over- 

shoot 

Max 

Fin 

Angle, 

Deg/G 

Max 

Fin 

Rate, 

Deg/s/G 

Gain 

Margin 

(dB) 

Phase 

Margin 

(Deg) 

Loop 

Gain 

Freq. 

(Hz) 

States 1 

Kcl 

0.385 

0.961 

-8.3 

-0.4 

0.928 

52.5 , 

-11,7.3 

± 33 

5.2 

2 ; 

Klqr 

0.238 

0.347 

-4.0 

4.2 

0.854 

15.0 

-6,+«> 

± 60 

3.7 

1 ? 

Koo 

0.213 

0.514 

-10.2 

0.0 

1-136 

61.4 

-5,6.3 

± 30 

4.8 

6 '[ 


0.329 

0.784 

-2.4 

0.0 

0.668 

8.0 

-8,10 

± 39 

3.4 

27 ' 

Krep. 

0.213 

0.784 

-2.4 

0.0 

0.668 

8.0 

-8,10 

± 39 . 

3-4 

22 

KooSF 

0.213 

0.514 

-9.3 

0.0 

1.053 

34.4 

-5,-l-co 

± 40 

24.0 

2 

KsPPC 

0.213 

0.514 . 

-9.3 

0.0 

1.053 

34.4 

-6,*fo® 

± 40 

4.7 

4 

Kppc 

0.213 

0.514 

-9.3 

0.0 

1.053 

34.4 

-6,-h» 

± 40 

52.1 

2 


* This column represents controllers designed using Classical, Linear Quadratic, 
DGKF Algorithm, p-Synthesis, Reduced Order p-Synthesis, Hoo State Feedback, Strictly 
Proper Projective Control and Proper Projective Control techniques, respectively . 
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2.8 Conclusions 

In this chapter parametrized longitudinal autopilot for missile has been designed 
using pole-zero assignment. The design technique employed Youla parametrization of 
stabilizing compensator which facilitated incorporation of poles and zeros through 
polynomials. The crux of the design is finding out the overall closed loop transfer 
functions which meet the design specifications . Direct search methods of nonlinear 
optimization have been found very suitable for selection of the closed loop transfer 
functions. For pitch rate loop ITAE optimal transfer function under the constraint of 
L2-norm on u(t) was chosen. For outer acceleration command tracking loop, objective 
function having the terms corresponding to tracking error, maximum overshoot, 
undershoot, gain and phase margins has been minimized to achieve the overall closed 
loop transfer function. Once these transfer functions are selected, the rest of the 
design is straight forward. 

The element of subjective judgement or requirement of hit and trial is much 
less as compared to other techniques. The design of parametric controller, in the 
frame work of algebraic theory, using optimization techniques is very potent area. 
The designer’s level of subjective judgement on the design reduces. The design is 
more transparent and the algorithms lead to controller from desired specifications in 
a systematic way. 

Disadvantages of the methodogy are (i) order of compensator is high (ii) use 
of conventional direct search methods for nonlinear optimization does not guarantee 
global optimum solution but only local optimum. But this is not a serious disadvantage 
in the sense that the solutions finally achieved have been found reasonably insensitive 
to the starting solutions. In fact even where the starting solutions were chosen well 
beyond dominant pole configuration, the final solution converged almost to the same 
value. Regarding high order of controllers, the model reduction techniques can be 
applied to reduce the order of the designed compensator. 



Chapter 3 


A Two Degrees of Freedom Parametric 
Longitudinal Autopilot Synthesis using 
Model Matching and Constrained 
Optimization with Simulated Annealing 

3.1 Introduction 

The objective of the control system design is to select the controller such that 
the closed loop system is stable and has the desired characteristics with respect to 
the reference command, disturbance, sensor noise, sensitivity and robust stability for 
parameter perturbations. In control systems with one degree of freedom configuration, 
it is not possible to cope up simultaneously with the two problems of : 

(i) attaining a desired system response with reference to command tracking 
and 

(ii) achieving feedback properties to handle plant parameter variations and/or 
disturbances acting on the plant [40]. 

It has been pointed out by Desoer et al [41] that properties with respect to 
the reference command and closed loop characteristics could be independently specified 
using a two degrees of freedom (TDOF) structure. 

In this chapter, longitudinal autopilot has been designed using parametrization 
of stabilizing controllers with two degrees of freedom configuration given by Hara 
[42]. This parametrization gives possible classes of stabilizing controllers with two 
degrees of freedom configuration using fractional representation approach. The possible 
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classes of controllers and realizable characteristics are provided by two free independent 
parameters of proper stable rational matrices. This two degrees of freedom configuration 
is for the case in which the measured variables are not necessarily coincident with 
the controlled outputs. The missile pitch dynamics form the similar system with one 
controlled output (pitch acceleration) and two measured variables (pitch rate and pitch 
acceleration). 

This parametrization of stabilizing controllers with two degrees of freedom 
configuration is superior to others given by Yoshikawa et al [43] and Vidyasagar 
[44] in the sense that the properties with respect to the reference command depend 
only on one parameter, Ari(.s), and the closed loop properties depend only on another 
one, K2(s). This facilitates independent shapings of time domain and frequency domain 
responses by proper selection of Ki{s) and K 2 (s), respectively. All the objectives of 
the autopilot design can be expressed as the hard bound constraints on time domain 
and frequency domain responses. Hard bound constraints on the time domain are in 
terms of the rise time, peak overshoot, undershoot, settling time, steady state error 
and on the control input. The main frequency domain hard bound constraints are 
related with disturbance rejection over the bandwidth of the feedback system and 
stability robustness to plant uncertainty and parameter variations. Stability robustness 
constraints are equivalent to the classical gain and phase margins. 

In this chapter, the time response shaping has been carried out using closed 
loop model maching concept [45,46]. This technique reduces the nonlinear optimization 
problem to a problem of linear simultaneous algebraic equations in terms of solution 
variables. 

The frequency domain response shaping has been carried out by casting the 
design problem as constrained optimization problem which has been solved using 
Simulated Annealing (SA) technique [52]. Constrained frequency domain optimization 
problem requires tuning of bounds on the performance functions during optimization 
process. This approach of frequency response shaping is superior to unconstrained 
Hoo minimization approaches [47-49] because the latter suffer from the following 
drawbacks : 
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p-- 

^ R(s) c,(s) —to ! 

"1 * 1 

L 1 

1 

1 

vjs) 1 

' ► M ^ j 

- r ; 


G(.s) 


Figure 3.1 : Feedback Control System with Two Degrees of Freedom 

Xis) N(s) + y(5) D(s) = I,n 

(ii) rank [N^(k) D^(X)] = m; V e Ct 

where Ct denotes the extended right half plane, ie. {5 : 5 > 0 } kJ {=« }. In 

this case, ND~ ^ is referred to as right coprime factorization of H. The left 
coprime factorization is defined analogously. 


3.2.2 Two Degrees of Freedom Configuration 


The two degrees of freedom configuration of the control system, given by Kara 
[42], is shown in the Fig.3.1. C(5) = [Ci(5) Ciis)] e ^ and 

R{s) = ol{s) 7 = 6 R§^^ denote the plant, the controller and the reference 

command generator, respectively. The signals r, e, u, ym and yc represent the 
reference, error, control input, measured output and controlled output, respectively. 
The matrix M is assumed of the form M = [Ip 0] and the total plant is defined by 


d - 


G{s) = 


MP{s) 

P(s) 


With the internal model for the reference command at the error channel, 

the servo system is shown in the Fig.3.2. Consider the augmented system of the 
total plant G(s) and the internal model //(j)(s) which is represented by [42], 
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Figure 3.1 ; Feedback Control System with Two Degrees of Freedom 
X{s) N{s) + Y{s) Dis) = 1,„ 

(ii) rank [N^(k) d'^{X)] = m; V X 6 

where Cj denotes the extended right half plane, ie. {5 : 5 > 0 } 1^ {<=0 |. In 

this case, ND~ ^ is referred to as right coprime factorization of H. The left 
coprime factorization is defined analogously. 


3.2.2 Two Degrees of Freedom Configuration 


The two degrees of freedom configuration of the control system, given by Kara 
[42], is shown in the Fig.3.1. C(^) = [Ci(5) Ciis)] e and 

R(s) = a(s) I = e RS^^ denote the plant, the controller and the reference 

command generator, respectively. The signals r, e, u, y,n and yc represent the 
reference, error, control input, measured output and controlled output, respectively. 
The matrix M is assumed of the form M = [Ip 0] and the total plant is defined by 


d - 


Gis) = 


MPis) 
P(s) • 


With the internal model l/^{s) for the reference command at the error channel, 
the servo system is shown in the Fig.3.2. Consider the augmented system of the 
total plant G(s) and the internal model I/^(s) which is represented by [42], 
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Figure 3.2 : Servo System with Two Degrees ot Freeiioiii 


d(s) 



O' 


m(s) 

Figure 3.3 ; Standard Feedback Configuration 
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Ge(s) 


MP{s) /(t)(5) 
P{s) 


(3.1) 


The C(^) = [ti(s) / (})(5) 6(^)] is a servo controller for G{s) if 


d 

t^{s) = [^iCs') ^2(-s)] stabilizes the augmented system Geis) under the assumptions 

of solvability condition (A2.5) or (A2.6) (given in Appendix 2), which assures no 
unstable pole-zero cancellation between G{s) and //(1)(5) [42], Hence the class of 
controllers t(s) that stabilize Ge(s), is investigated. 


3.2.3 The Class of Stabilizing Controllers 

In the following, the class of stabilizing controllers is given by the Lemma 3.1 
which is taken from [42]. Other details pertaining to the parametrization of controllers 
for the present TDOF configuration, which are based on [42], are relegated to 
Appendix 2. 

Lemma 3.1 

Suppose that solvability condition (A2.5) or (A2.6) holds. The controller 


^(s) = [&i(^) ^2(‘S')] which stabilizes Ge(s) is given by 

&i(i) = (Xi + i>Ki) , (3.2) 

tiis) = (V-K 2 ^-^ a-Ki X 2 ]+ K 2 D) (3.3) 

for any Ki(s) E and K 2 (s) E R-^^. 

Furthermore, the attainable transfer characteristics are expressed as 

Her{s) = (3.4) 

Hyr{s) = ^ {Xi+i^Kx) ■ (3.5) 

Hur{s) =XD\ + Dbm (3.6) 

Si{s) = D b {Y - KiNy (3.7) 

<?0'\-bK2) 

Y2-N2bK2 


Sq{s) = 


(3.8) 
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where Her(s), Hyr{s) and Hur(.s) are the transfer function matrices from r to e, 
r to yc and r to u, respectively, and Si{s) and So{s) represent the relative sensitivity 
matrices at the input (u) and the output (y,„) channels of the plant, P(s), respectively. 


Remark 3.1 

Setting <^(s) = 1 in Lemma 3.1, we have the results for TDOF regulator systems. 

Remark 3.2 

It can be noted from equations (3.4) - (3.8) that the properties with respect to 
the reference command and sensitivity characteristics depend only on Ki(s) and K2(s), 
respectively and these two types of characteristics can be specified independently in 
the control system design. 


3.3 Design For Pitch Acceleration Command Tracking - 
using Closed Loop Model Matching Concept 

3.3.1 General 

Closed loop system in pitch channel is required. to track reference acceleration 
command. The typical time response specifications for a command tracking system 
are rise time, peak overshoot, undershoot, settling time and steady state error. Moreover, 
the command input signals can drive the missile plant outside the linearization region 
which may lead to performance deterioration or even instability. So design specifications 
should include some measure of maximum actuator deflection and maximum actuator 
rate. 

In TDOF parametrization used here, the reference command tracking depends 
on parameter isri(5). Thus the design objective is to select proper ^^ 1 ( 5 ) which will 
give the desired tracking performance. The closed loop model matching concept 
[45,46] has been used here for determining ifi(s). 
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3.3.2 Problem Formulation 


In this Section, the reference command tracking problem is formulated as closed 
loop model matching problem. With internal model present in the Fig. 3.2 for the 
reference command in two degrees of freedom configuration, the equation (3.5) is 
rewritten as 

Hyr(s) = k iX\+m) . (3.9) 

The cost function in closed loop model matching problem is defined as a mean 
square error function [45,46] 


J = 


^2 

1 

J 


j^rO'“)-HyK/(o) 


2 

do) 


(3.10) 


Here, H^r(j O) ) is the desired closed loop transfer function, (coi , C 02 ) is the 
frequency interval of interest over which the closed loop model matching is to be 
achieved. 


Let 


Hyrf/CO) = 


Z2 


HyrO'®) 


Hi _ Ri+Ry 
Hi Si + Sij 


Ni(j(o) 

XiO'O)) 

Kim 


Cl + Cij 

+ gv 

di + dy 
gl + gzi 
Ai + Ay 
Bi + By 


i+s” 

1 = 1 


(3.11a) 

(3.11b) 

(3.12) 

(3.13) 

(3.14) 


Here, the aim is to solve for unknowns aCs and bi's. One approach may be 
to set VJa,b = 0 and solve a set of nonlinear simultaneous equations. In this approach, 
rather than minimising equation (3.10) ie. the quantity 
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^2 


I (^HiZ2-H2Zi^ /H2Z2 p 4co, 


coi 


^2 


the quantity 1 H 1 Z 2 - H 2 Z 1 

coi 


d(£) is minimized. If 


a solution exists, the 


minimum value of this quantity will always be zero. Therefore, under optimal 
conditions H 1 Z 2 = H 2 Z 1 , provided that polynomials on both sides of this 
equation are of the same order. 

Upon substitution of equations (3.9) and (3.11) - (3.14), and multiplying the 
resultant integrand by its common denominator, the following equation is obtained. 


where. 


r - 


j = 


B\C- B 2 D - PA\ + QA 2 


+ {B\D\B2C-PA2-QA\ 



4© 


COI 

(3.15) 

2 4 

Al = aO“a2CD +04© 

(3.16) 

3 5 

A 2 = 01 © -03© + 05(0 -.... 

(3.17) 

Bi = bo-b2(i>^ + b4(i)^- ... 

(3.18) 

B 2 = bid) - b3(£i^ + bs(£i^ - ... 

(3.19) 

P = (SiCi- S 2 C 2 ) 51 - (5iC2 + 52 C 1 ) g 2 

(3.20) 

Q = (SiC2 + 52Ci)^ 1 + (5lCl -52C2)g2 

(3.21) 

C = E-G 

(3.22) 

D = F-H 

(3.23) 

E = Rlgl 82 

(3.24) 

F = 2R\ g\ 82 + R2(81- gh 

(3.25) 

G = (5i Cl - S 2 Ci) d\ - {S\C 2 + 52Ci) d 2 

(3.26) 

H = (Si C 2 + S 2 Cl) 4i + (SiCi - 52 C 2 ) d 2 

(3.27) 


On setting 

V Ja.b = 0, 

a set of linear algebraic equations is obtained. These equations can be written 
in the matrix form as 
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JiT = S (3.28) 

whe re. 


Uo 

0 

- Ui 

0 


Wi 

V2 

- W 3 

-V 4 .. 

0 

U 2 

0 

- c/4 


- V2 

W3 

^4 

- 1^5 •• 

- Ul. 

0 

c /4 

0 


- W3 

- V4 

W5 

V 6 .. 

0 


0 

t/6 

......... ... 

V4 

- W5 

-V6 

W7 .. 


Z2 0 -Z4 0 ... 

0 Z4 0 -Ze ... 

-Z 4 0 Z 6 0 ... 

0 -Z6 0 Zg ... 


(3.29) 

" Vo ■ 

S W\ 

; - V2 


: 0 
: Z2 
- 0 
i 24 


L • - (3.30) 

where , 

r 

Vi = 1 0)' (P^ + Q^) dco 


«0 

ai 

a 2 

as 



^2 

b4 



-V2 

- W 3 

74 


^3 

-74 

- W 5 

-T3 

74 

Ws 

-76 

- V 4 - 

~W5 

76 

W 7 


COl 


(3.31) 
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Wi = 


^2 

I co' {PD - QQ doi 

) 

COl 


Vi = 


^2 

I (s)‘ {PC + QD) da 
(01 


( 3 . 32 ) 


( 3 . 33 ) 




Zi = 


co' {C^ + D-) da 


COl 


( 3 . 34 ) 


Ari( 5 ) can be obtained for the closed loop model matching problem by solving 
the above set of equations. The above formulation considers only acceleration command 
tracking problem. A large control signal at actuator may have large actuator deflection 
and if the rate of rise of control signal is large, it may produce large actuator rate 
which may drive the missile plant out of linearity zone and even may lead to 
instability. 

This problem can be eliminated by introducing in the cost function the 


^2 

term | | Hur P do, which may be interpreted as square of L2-norm of the control 

(01 

signal due to a unit impulse within the frequency range (coi , (02)- The modified cost 
function can now be written as. 


^2 

/ = I I Hyr-Hfr 

J 


^2 


c/CO + 1 -^ur 

J 


“1 “1 ( 3 . 35 ) 

Upon substitution and after some manipulations on earlier lines, cost function reduces 


to 



A TWO DEGREE ... ANNEALING 


44 


^2 


/= I 
J 




PA\+QA2 + BiC-B2D'^ +^-PA2-QAi+BiD + B2C'' 


COl 


+ k 


2 


(^d\B\-d2B2 + g\A\-g2A2'j +(^d\B2 + d2B\+g\A2 + g2A\'^ Jl-ico 


(3.36) 


On setting V fa, £? = 0, the equations in the same form ie XY = Z are obtained 
with modified Uu Wi, Vi and Z/ ‘s. In this case, these are given by 


^2 

Ui,m = I + + 


4co 


COl 


(3.37) 


^2 

Wi,m= I co‘[(/>D-!2C)+/:(^-5i42 + ^24i j] 

COl 


4co 


(3.38) 


^2 


Vim = CO 


‘j^ ( PC+ QD) +^^-^l4l -^ 2 ^ 2 ^ 4co 


COl 


(3.39) 


^2 

Z/.m = I co‘ + + + 

COl 


d(i) 


(3.40) 


X, Y and Z will have the same form as given in equations (3.29) and (3.30) 
with the difference that Ui, m, Wi, m, Vi, m and Z/, m, will replace the respective terms 
Ui, Wi, Vi, and Z/ in equations (3.29) and (3.30). 

By varying the value of k, a tradeoff is established between reference acceleration 
command tracking and control signal (u), and its rate («). 


3.3.3 Algorithm 

For implementing the above approach the following algorithm is developed for 
selecting Xi( 5 ) to achieve the desired tracking response: 
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Step 1 

Select a second order transfer function that gives the desired time response 
characteristics. Take this as the desired transfer function. 

Step 2 

Determine the plant and its Bezout coprime factors. 

Step 3 

Set = ^k\ni = 0 and k = 0 {6 denotes the degree of polynomial) 

Step 4 

Check, if 5 (Hi Z 2 ) = 5 (H 2 Zi). If yes, go to Step 5; otherwise increment 
8kin and ^ (Hi Z 2 ) = 5 (H 2 Zi) is satisfied. 

Step 5 

Check if Ki(s) is stable. If yes, go to Step 6; otherwise increment dku and 
and go to Step 4. 

Step 6 

Evaluate tracking response. If satisfactory, go to Step 7; otherwise increment 
and and go to Step 4. 

Step 7 

Evaluate actuator deflection and actuator rate. If satisfactory, stop search for 
otherwise go to Step 8. 

Step 8 

Take k ^ 0 , and solve for Ki(s) using equation (3.28) with C//, nil Vi, /M, 
Wi,fn and Zi,m, as obtained from equations (3.37)-(3.40). 

Step 9 

Evaluate tracking response. If satisfactory, stop search for .^ 1 ( 5 ); otherwise modify 
k and repeat Step 8 until satisfactory tradeoff is obtained between tracking response and,l 
actuator deflection and actuator rate. If satisfactory trade off is not obtaine4 increment 
bkin and him and go to Step 4. 
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The problem is posed as constrained optimization problem and is solved using 
Simulated Annealing (SA) technique. SA technique is a global* optimization technique 
which finds out global optimum in the presence of many local optima. 


3.4.2 Frequency Domain Performance Specifications 
(i) Disturbance Rejection 

Disturbance and measurement noise rejection can be obtained by keeping 
\Si (/oo)| small at lower frequencies and \Ti (/ol))| small over higher frequency range, 
respectively. So good disturbance rejection over the bandwidth of the feedback system 
can be achieved by keeping the |5i (/m)| small over the bandwidth (0, (Oc)- 

Hence we define the following performance function to capture the above 
mentioned constraint. 

(3.41) 

where &/(co) is a continuous bound function. 

Therefore, the disturbance rejection performance requirements can be formulated 
in terms of the following inequality 

^i(k2) < 0 (3.42) 

A simple choice for the bound function would be to set [49] ; 

0 < 6/(cd) = bfi « I , if CO < (Oc (3.43) 

£)/(co) - bf 2 > I , if (0 > cOc (3.44) 

It follows from the extension of Bode’s integral theorem to multivariable systems 
[49,54,55], that if the sensitivity is made lower at some frequencies, the penalty is 
a higher peak value of the sensitivity at some higher frequency. So &/(co) must 
exceed over some frequency range beyond cOc 
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Some time, the time response specifications of the selected transfer function 
may not be achievable even after a large number of iterations. In such case, modify 
the desired transfer function and repeat all the steps. 


3.4 Design For Shaping Frequency Domain Responses 
3.4.1 General 

Feedback in a closed loop system is required to achieve various desirable 
properties, such as stability, disturbance attenuation and low sensitivity to changes 
in the plant. Since all these properties depend on the shape of various responses in 
the feedback loop, the simultaneous loop shaping of several frequency domain responses 
is essential. 

In the present approach of frequency response shapings, constrained optimization 
technique has been used. The various frequency domain hard bound constraints on 
the performance functions are defined and the problem is posed as constrained 
optimization problem. These bounds are tuned during optimization to achieve satisfactory 
performance. The main frequency domain hard bounds are related with disturbance 
rejection over the bandwidth of the feedback system and stability robustness to plant 
uncertainties and parameter variations. If no hidden modes are present, then the closed 
loop stability can be inferred from the stability of either the output or input sensitivity 
functions [53]. In missile pitch dynamics, no hidden modes are assumed to be present. 
So the loop shaping has been done with the loop broken at the input of the plant. 
Since the pitch dynamics form a Single Input Multi Output (SIMO) system, with 
one input and two outputs, the loop transfer function, when the loop is broken at 
the input of the plant, turns out to be scalar transfer function. Robust stability 
constraint has been imposed in terms of M-circle constraint [51], which provides 
robustness against parameter variations and sector bounded time varying nonlinearities. 
M-circle constraint guarantees given classical gain and phase margins. Apart from 
the above frequency domain constraints, other constraints are incoporated to capture 
the desired sensitivity .and complementary sensitivity functions responses. 
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Loop gain 
Margin Constraint 




Real part of Loop Gain- ► 

Figure 3.4 : Diagram showing M-Circle con.straini 

(ii) M-circle Constraint For Robust Stability 


M-circle constraint provides an useful way of specifying stability margins [51], 
The loop margin, M, is simply the minimum distance in the complex plane between 
the Nyquist plot of the loop gain and the critical point (-1 +j0) : 


M 


min 

CO 


1 + L/(/co) 


(3.45) 


Here, Li (/co) is the loop transfer function with the loop broken at the input 
of the plant. The M-circle constraint requires the Nyquist plot of L/ (/(U) to maintain 
at least a distance of M (the M-circle radius) from the critical point (-1+jo), as 
shown in Fig.3.4. This will provide guaranteed gain and phase margins. 

Therefore, a classification of minimum loop margin, 

M > A/min (3.46) 

can be interpreted as robustness specification [51], as shown in Fig. 3.5. The 
perturbed loop gain is L{s) = Li{s) + A(s), where A(s) is a stable transfer function. 
Such a perturbation is called additive loop gain perturbation. 
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Figure 3.5 ; The Perturbed Plant for Loop Margin Constraint 
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Rgure 3.6 . System with a Time Varying Nonlinearity 
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The closed-loop system is stable for all stable A such that || A ||oo < . 

•^min 

In Other words, meeting a minimum loop margin specification is equivalent to 
guaranteeing stability of the closed loop system despite additive perturbations of the 
nominal loop gains by any stable transfer function with magnitude less than 1/M for 
all frequencies. 

The loop margin, M, is related with the sensitivity function as follows [51] : 


M/Z 1 l+L/(/'co) 1 

(3.47) 

1 


Max 

(0 1 + L/go) 

,(3.48) 

1 


||5;(^)ll~ 

(3.49) 


Thus small loop margin corresponds to peaking of the sensitivity function at 
some frequency. The equation can also be expressed as 


II Si W II ~ 



(3.50) 


using the inequality given by equation (3.46). 


Another interpretation of M-circle constraint can be based on circle criterion 
which provides sufficient condition for closed loop stablility despite certain nonlinear 
plant perturbations [51]. The perturbations are time varying memoryless nonlinearities 
in the actuator which can be expressed as w = t) (Fig.3.6). The time varying 
memoryless nonlinearity f(.) satisfies a sector condition : roughly speaking, _/(.) 
multiplies u by at least a and at most b, i.e. for all t, 

ax^ < xf(x,t) < bx, (3.51) 

where we assume 0 < a < I < b. This is shown graphically in Fig.3.7. 

The circle theorem states that if the nominal system is closed loop stable and 
the Nyquist plot of the nominal loop gain remains outside the circle symmetric with 
respect to the real axis and passing through the points 1/b and 1/a, then the perturbed 
closed loop system of the Fig.3.6 will be stable. 

The circle theorem yields another interpretation of (3.46) : constraints expressed 
by equations (3.46) and therefore (3.49) hold if and only if the circle criterion is 
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Figure 3.8 : Topical gain Characteristics of sensitivity and 
complementary sensitivity functions 

extends over a much wider frequency range. Therefore, \Si0(D)\ is kept small at 
low frequencies (for 0 < oa < CDo, say) and |r;(/0))| is kept small at high 
frequencies (for co > (Oi , say). 

Though performance functions (t)i(ii:2) and already defined, both put 

constraints on \Si{s)\ in various frequency ranges, a few additional performance 
functions are required to cover more frequency points/ranges and capture the typical 
shapes of |5'i(/a))| and |7i(/(jo)| as shown in Fig.3.8. 

In a very low frequency interval 0 < O) < coo, (where coo is much lesser than 


(He) we require that : 



\Si(j(H)\ < pi. 

for 0 < (0 < too 

(3.54) 

and we have 



|%w)| = 1, 

for 0 < CO < coo 

(3.55) 


Here pi (« 1) is the bound on |.5/ (/(o)| in the frequency range 0 < co < (Oq. 

We define the performance functions ^a{K' 2 ) corresponding to equation 

(3.54) and (3.55) as : 


(k2) = 


Sup 

CO e [0,coo] 


{\Sim\-piJ 


A Sup 

~ CO e [0,coo] 


(3.56) 


<J)4 (Kl) 




(3.57) 
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Figure 3.7 : An Example of a Sector [a.b] Nonlinearity f(x,t) 
satisfied for nonlinearities in sector [1/(1 +M), 1/(1 -M)]. Thus constraint (3.49) 
can be viewed as guaranteeing stability despite time-varying memoryless nonlinear 
perturbations in the loop. 


Thus, the M-circIe constraint performance specification can be defined as 

1 


[ 0 ,-) 


5, (/CO) 


Atrain 


Therefore, for meeting the M-circle constraint we require that 
<^2{K2) < 0 


(3.52) 

(3.53) 


(iii) Input Sensitivity And Complementary Sensitivity 
Shapings 

It has been shown in [48] that many essential design requirements can be 
achieved by shaping of sensitivity and complementary sensitivity functions. The 
unavoidable trade off between these two is given by equation (A3. 14) of Appendix 
3. In practice, the conflict between keeping T(s) small and keeping S(s) small is 
resolved by making one small at some frequencies and the other small at other 
frequencies. Usually, the spectrum of reference signals and disturbance signals are 
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Here pi is constant and provides bound on the \Ti (/a))| in the frequency 
interval (0, (Uo) 

The bandwidth of the feedback loop ,0)^, is defined as the lowest frequency 
such that 

|7'/(/'®t)| = mi/^ = 0.707 (3.58) 

Performance function corresponding to the value of Ti (/co) at bandwidth is 
defined as 

A 

<^5(i<:2) = (|r(/coi)|- 0.707) (3.59) 

The loop bandwidth, (Db, is usually very close to ’0 dB crossover frequency’ 
(Oc of the loop transfer function Li(s) (= C(s) G{s)). For acceptable designs, loop 
bandwidth can be estimated in terms of crossover frequency by the relation[32]: 

(Dc < o)i, < 2 (Oc. (3.60) 

Therefore, we define, correspondingly, a performance function 

<^6 (Ki) = {\L(Joic)\-l} (3.61) 

Permissible peak of |5/ (j co)| has been specified in performance function 
(Ki) and also in the M-circle constraint ^i (Ki). Permissible peak of the 
\Ti (/(o)| is specified by another performance function : 

(3.62) 

where, pt is the bound on \Ti (/(o)| in the frequency interval (O £ [0, oo/j]. 

Another performance function is required to capture the shape of \Ti (j oS) \ 
beyond the frequencies greater than (Db . It is required that 

1 Ti(j(D) \ > ph , for (0 > m. (3.63) 

The corresponding performance function is defined as 

rnmi-pu 

f 

Therefore, to capture the shapes of \Si (j (o)| and \Ti (/(o)| at various frequencies, 
we require. 


93 (/!: 2 ) < 0 


(3.64) 
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b4 (Ki) < 0 

(3.65) 

b5 m = 0 

(3.66) 

(|)6 m = 0 

(3.67) 

b7 m < 0 

(3.68) 

(i>8 (^ 2 ) < 0 

(3.69) 


3.4.3 Objective Function 

All the design requirements in the design for the frequency domain loop 
shaping are expressed in the form of either equalities or inequalities . Weighted sum 
objective function with external penalty is used in the following form : 

A 

(K 2 ) = /S, < <Pi (K 2 ) > ^ + h [ <Pj {Kl) f } 

where cp/ (Kj) and cpy (K 2 ) are the performance functionals with inequality and 
equality constraints, respectively. In the first term on the RHS of the above perfomance 
function , < a > = a , when inequality is violated; it is zero, otherwise. Second term 
on the RHS is used to handle equality constraints. Since the feasible point always 
satisfies (pj (K 2 ) — 0 , any infeasible point is penalized by an amount proportionate 
to the square of the constraint violation . Xi and Xj are the non-negative scalar weight 
factors. 

During the optimization process, the Xi and Xj, and bounds are tuned to achieve 
the optimal solution. Here, <^ 2 (K 2 ), ^3(i^2), MKl), and iK 2 ) are the 

performance functionals with inequality whereas ^5 (K 2 ) and <1)6 iK 2 ) are the performance 
functionals with equality. 

3.4.4 Simulated Annealing Technique 

Simulated Annealing (SA) is a global optimization technique that distinguishes 
between different local minima. Starting from an initial point, the algorithm takes a 
step and the function is evaluated. When minimising a function any downhill step 
is accepted and the process repeats from this new point. An uphill step may be 
accepted. Thus, it can escape from local minina. This uphill decision is made by the 
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Metropolis Criterion [52,58]. As the optimization process proceeds, the length of the 
steps declines and the algorithm closes in on the global minimum. The SA algorithm 
[52] with some modifications as in [63], has been used for optimization, details of 
which are given in Appendix 4. The SA algorithm used here is applicable to the 
functions which do not need to be smooth or even continuous in their domain. 


3.5 Example 

Consider the same missile pitch dynamics model as taken in Example 2.6 of 
Chapter 2. 

Typical design objectives for the autopilot can be stated in terms of the time 
and frequency response specifications as follows: 

Design an autopilot to track the commanded acceleration maneuvers with a 
(63%) rise time less than 0.3 secs, overshoot less than 5%, undershoot less than 
2.5%, (95%) settling time less than 1 sec. The autopilot must provide the gain 
margins of (-6, +10) dB and phase margins of ± 40 degrees. The autopilot bandwidth 
is to be restricted by the high frequency unmodeled dynamics constraints (e.g. 
unmodeled flexible modes and actuator nonlinearities). 

The present design approach facilitates independent time response and frequency 
response shapings. Now time response and frequency response designs are taken up 
separately to determine iiri(s) and K 2 (s), respectively. 


3.5.1 Design For Acceleration Command Tracking 

For the closed loop model matching, the desired closed loop transfer function 
is chosen of the second order given by, 

^ ^ 55 

(^^+ 10 . 5 ^ + 55 )’ ( 3 . 71 ) 

which gives the desired time response. 

The left and right coprime factors and the corresponding elements of the Bezout 
identity, are obtained by the algorithm of Nett et al [59] which requires to solve 
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only two pole placement problems. These doubly coprime factors are the solution 

of the following equation 

Yis) Xisf'Dis) -X(s)~ '/o' 

-Nis) D{s) N{s) ?(s) _0 /_ (3.72) 

These two pole placement problems have been solved with arbitrary pole 
locations at ^7.8 i 90. 4i, —10.713, —1.563. 

The stable coprime factors relevant to the design of time and frequency responses 
are given below : 

G = /+ 147.885^ + 14450.0/+ 1.5902e5j + 2.1381e5 (3.73) 

’ (3.74) 

where 

x\ = 8.8098e - 3 1.2195^+ 1.1852e2^ + 7.745e2 


where 


where 


m = -2.4254e6/+ 1.6644e9 



4 = / + 136.9j^ + 128985^ + 101895 - 609220 


(3.75) 


(3.76) 



where 


= s‘^+l.58E5e2s'^+L6l23e4s^ + 3.46Se5s + 2.l4S5e6 
N\ n\ 

=r - ’ 

' N 2 «2 

where ni and n2 are given by : 

m = -2.4254e65^ + 5.6313e65 + 5.6372e8 
n2 = - 1.4693e65- 1.5737e7 


Du Dn" , 3n 3i2 

D = - — 

^ ~ G ’ 

D 2 \ D 22 321 322 

where 3ii, 3i2, 321 and ^22 are given by : 


(3.77) 
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311 = /+ l,415e2 5^ + 1.3328e4j-^ + 4.9972e4j + 6.6777 e5 

312 = 2.9054e3i^ + 3.8518e5 + 3.0566^7 5 + 2.7443e7 

321 = 3.9573e-2 5^ + 7.35865^ + 6.1302e2 j + 5.5711e3 

322 = / + 1.4328e2 + 1.4107e4 ? + 1.35S4e5 j + 3.3886e4 

For regulator systems, (p (s')=l, and for the class of step reference input signals, 

the solvability condition given by equation (A2.5) is satisfied for all cp ( 5 ) = 
if no root of ni is at the origin and G does not have any root at s = - a. 

Now decomposing the factors into real and imaginary parts as defined in 
equations (3.11) - (3.13), we get 


Cl = 2.4254e6co^+ 1.6644e9 

(3.78) 

11 

0 

0 

(3.79) 

^1 = 14450 co^ + 2. 1381 e5 

(3.80) 

g2 = 1.5902e5(O- 147.88 0)^ 

(3.81) 

3i = 7.7475e2- 1.219 CD^ 

(3.82) 

32 = 1.1852e2(D-8.8098e-3co^ 

(3.83) 

Ri = 55.0 

(3.84) 

0 

II 

(3.85) 

Si = 55 - co^ 

(3.86) 

52 = 10.5 CO 

(3.87) 


These equations are substituted in equations (3.20) - (3.27) to evaluate the 
integrals given by equations (3.31) - (3.34). Ui,Vi,Wi and Zj’s are found out by 
evaluation of integrals given by equations (3.31)- (3.34) in the frequency interval 
(03i , C 02 ). The solution vector, Y, in equation (3.28) is 

Y = [ao , ai , a2 , , bo , bi , b2 , These a/’s and bi’s are the coefficients of 

numerator arid denominator of Ki(s), respectively. The matrix X and vector Z in 
equation (3.28) have entries made of Ui , Vi , Wi and Zfs. The order of the system 
of equations depends on the number of at's and bfs to be found out for the selection 
of Ki{s). The number of a/’s and bfs depends on the degrees of the numerator 
{hiJ and denominator (6yt,„) of Ki{s). 
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Since ^^21 (5’) and K 22 {s) are proper stable rational functions, the numerator and 
denominator will have the same degree i.e. m=n. 

The gain and phase margin constraints are imposed by the M-circle constraint 
given by equation (3.50). Mmin = 0.7 guarantees gain margins of (-4.6, +10.4) dB 
and the phase margins of ± 41 degrees, approximately. Initially the design is aimed 
at achieving C0c=30 rad/sec and o)z>=50 rad/sec. 

For frequency shaping design the objective function given by equation (3.70) 
is minimized to get solution vector, {a'l ,b'i ,c'i} using Simulated Annealing 

algorithm given in Appendix 4. Scalar weight factors have been used in the 

objective function. Initially weight factors, X/, z=l, ,8 , were taken unity. The 

following initial values of the bounds were taken: 
bfi = 0.4, at (0 < cOf 

bf2 = 1.4 

pi = 0.2, at coo = 5 rad/s 
P2 = 1.0 
pt = 1.5 

The following parameters were used for running Simulated Annealing algorithm 
given in Appendix 4: 

Ns = 20 
Nt = 60 

Ci - 2 , i = 1, .... ,n 
Ne = 4.0 
IT = 0.85 

The SA algorithm was run for various initial guesses of the design variables 
but irrespective of initial guess solution converged. The starting temperature of l.SelO 
gave optimal solution. The implementation of SA algorithm given in Appendix 4 was 
run. Only one run of the SA algorithm was found necessary to get converged solution. 
However, the convergence of the algorithm dependend on the initial temperature 
selected. The above temperature value was selected using the criteria described in 
[52]. This value of temperature was arrived at after mnning many trial runs. Once 
the temperature parameter was selected, the program was run for the design variables. 
After each run, the performance of the design was evaluated. Accordingly, the 
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non-negative scalar weights and bounds were modified to get satisfactory design 
specifications. The following set of scalar weight factors and bounds gave satisfactory 
performance : 

Scalar weight factors : 

A.I = 42.0 
Xi = 1.0 
Xs = 11.0 
U = 300.0 
Xs = 350.0 
Xe = 30.0 
Xt = 50.0 
Xg = 50.0 


and bounds : 

bfi = 0.55 
bf2 = 1.3 
p\ = 0.4 
P2 = 1.0 


Pt = 1.4 . 

Ph = 0.7, at CO/, = coi, 

And the values of (Db and o)c were to be adjusted to 42 and 25 rad/sec 
respectively. This gives the following values of iif2l(5') and K 2 i{s) 

-1.11478e-2 + l.llllAZleS s + 2.974556e6 

(3.93) 


K2\{s) = 


5^ + 3005 + 22500.0 


K22{,s) = 


-4.584468e5 5^^ + 1.064291e6 5 + 1.06549519e8 
5^ + 300 5 + 22500.0 


(3.94) 


s 

With ({)(^) = and on substituting the values of Ki{s), K2\{s) and 

K22is), alongwith the values of X,Y,N,X2 and D the Ci( 5 ) and 
C 2 ( 5 ) = [ C21 (5) C22 (.y) ] , are given by 
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(5.32S6e-5 5^ + 2.4079^-2 + 4.4272 s'’ + 4.43376e2 / + 2. 1544^4 5^ +2.S027e5 / + 

1.5815e6j^+ 1.021e7 j- + 4.1629e7 j + 4.467Se7) 

' ” (7.i342e-2 / + 3.367el / + 4.658e3 s’’ + 5.4447e5 5^ + 2.0602e7 + 2.5557e8 / + ^ 

1.3692e9 + 3.939e9 + 5.9S76e9 ^ + 3.6226e9) 


(3.95) 


(-8.4859e-4/ + 1.8519e4^'^ + 2.9732e6 / + 3.0444e8 / + 6.6739e9 
7.216el0? + 5.4172eil + 2.3444el2 j + 2.5799gl2) 

■’ ~ (7.1342tJ-2/ + 3.3487el/ + 4.5722e3 / + 5.3275e5 / + 1.9236e7 5'* + 

2.0627e8i^ + 8.4055e8 + 1.7847e9i + 1.4134e9) 


(-3.2707e4/ - 4.954Se6j’ - 4.3437e8 5^ + 1.8414e9 / + 8.6124ell + 

^ , 1.6617el3 5^ + 1.0232el4^2 ^ 1.7712ei4i + 8.524el3) 

ClliS) = Q 1 A ^ 

(7.1342e-2/ + 3.3487el / + 4.5722e3/ + 5327 5e5 + \.9236e7 + 

2.0627e8? + 8.4055e8j^ + 1.7S47e9j + 1.4134e9) 


(3.96) 


(3.97) 


The following design specifications are achieved using the present frequency 
response shaping approach : 

Gain Margins = (-12.316, +14.8) dB 

Phase Margins = ± 82.2 deg 

Gain Crossover Freq. = 5.7932 Hz. 

The frequency response shapings are shown in the Fig.3.11 for Nyquist plot of 
Lj O’co), in the Fig. 3.12 for Bode Magnitude plot of Liija) to show loop shaping 
near the gain crossover frequency, and in the Fig.3.13 for gain characteristics of 
5j(/(o) and Tf (/co). 


3.5.3 Results and Discussion 

Achieved time domain performance specifications are shown in the Table 3.1 
and the frequency domain specifications are given in Section 3.5.2. The specifications 
achieved using present design methodology are comparable with the other designs 
reported in Table 2.2 of Chapter 2, which uses the same missile dynamics. Since in 
the present design approach, the time domain and frequency domain specifications 
have been achieved independently, any set of design specifications listed in Table 
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Figure 3.11 : NyquistPIot 
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Figure 3. 12 : Bode Magnitude Plot with the Loop Broken at the 

Input of the Plant 
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Figure 3.13 : Sensitivity and Complementary Sensitivity Functions at 

the Plant Input 


parameter k in the acceleration command tracking design, and bounds and scalar 
weight factors in the frequency response shaping. 


3.6 Conclusions 

The design approach developed in this chapter provides a systematic procedure 
for the acceleration command tracking design and frequency response design. Since 
the time and frequency response designs are independent to each other, wide variety 
of design specifications can be achieved through present design methodology. Simulated 
Annealing (SA) technique of optimization has proved to be very suitable for solving 
the constrained optimization problem. 


Chapter 4 


Synthesis of Decoupled Roll- Yaw 
(Lateral) Autopilot using Constrained . 
Optimization with Simulated Annealing 

4.1 Introduction 

Due to highly asymmetric BTT missile configuration, strong roll-yaw coupling 
is present. A BTT missile is first rolled to desired orientation and then required 
magnitude of force is achieved by adjusting the angle of attack ( a ) . If some 
sideslip is present, then in conjunction with the angle of attack large roll moments 
will be generated which are undesirable. So ideally, the sideslip should be zero. The 
lateral (roll-yaw) autopilot requires the closed loop system to command roll rate 
(stability axis) while keeping the sideslip small. 

In MEMO systems where input-output coupling is present, the decoupling 
methods [20, 33, 44] facilitate individual channel design. After decoupling, diagonal 
closed loop transfer function matrix is obtained.- 

As discussed in Chapter 2, in a control system the pole dominates transient 
response as well as stability [20-23] and the zero of a system plays major role in 
the interaction between the system and its external environment [24-26]. So wide 
variety of specifications can be achieved by assigning poles and zeros. After the 
input-output decoupling of the system, the poles and zeros can be suitably assigned 
to achieve the desired performance while maintaining stability of the closed loop 
system. 
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The 7/co-norm based constraints on the closed loop transfer function and sensitivity 
function have been found useful for specifying the robust stability against parameter 
variations and time varying nonlinearities [51]. M-circle concept [51], which is related 
with the sensitivity function, is convenient for imposing robust stability constraint 
against the plant parameter variations and the sector bound time varying nonlinearities. 

The roll-yaw channel of the BTT asymmetric missile considered here use roll 
and yaw rates for feedback. In its maneuver the missile is required to roll and then 
develop an angle of attack to intercept the target. So the roll rate is commanded. 
The roll rate commanded for the model under study is the stability axis roll rate 
p^. At the same time the lateral autopilot is required to minimise the sideslip ((3). 

The roll-yaw autopilot designed here is based on the decoupling approach 
studied in [60, 61]. The choice of the sensitivity function and hence the closed loop 
transfer function is made on the basis of time domain and frequency domain 
performance specifications e.g. tracking performance, sideslip, actuator effort, robust 
stability against parameter variations, uncertainties and time varying nonlinearities, 
disturbance rejection etc. The H^-novm based constraints are imposed and the problem 

is formulated as nonlinear optimization problem with constraints. This optimization 
problem aims at finding the coefficients of the polynomials involved in sensitivity 
and closed loop transfer functions which in turn will give parametrization of the 
controller. The problem so formulated is solved using Simulated Annealing (SA) 
technique [52]. The SA algorithm used here is found to give global optimal or near 
global optimal solution . Other details of SA technique are given in Chapter 3 and 
Appendix 4. 

The stmcture of the chapter is as follows : in section 4.2, the roll-yaw missile 
dynamics and closed loop design goals are described ; the control law based on 
decoupling and pole-zero assignment is discussed in Section 4.3. In Section 4.4, the 
time domain and frequency domain performance specifications are described and the 
design problem is formulated as optimization problem in Section 4.5. The design 
methodology and design algorithm are developed in Sections 4.6 and 4.7, respectively. 
The design algorithm is implemented in Section 4.8 and the Chapter is concluded 
in Section 4.9. 
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4.2 Roll- Yaw Missile Dynamics And Closed Loop Design Goals 

The roll-yaw missile dynamics involves the dynamical equations in sideslip (|3), 
roll rate (p), yaw rate (r), and actuator deflections. The linearized roll-yaw dynamics 
for the BTT missile, taken from [62], which is same as in [7], is described in 


Appendix 1. The governing equations are : 

P = psin(a)-rcos(a) + yp-P + ysa-5a + y5r-5r (4.1) 

p — Lp'P + Z.5 a • 5a + r • 5r (4.2) 

r == N^-^+ N&a-5a + N&r-5r (4.3) 

5a = - 2 ^ CO Sa - CO^ (5a - 5ac) (4.4) 

si- = - 2 ^ CO 5r - CO^ (5r - 5rc) (4.5) 


The roll rate (p) and the yaw rate (r) are measured outputs and have been 
used for feedback. In the state -space representation, the missile dynamics, in roll-yaw. 


can be expressed as : 

X = AX+BU (4.6) 

Y=CX (4.7) 

where, 

X = [P,p,/-,(!),v,5a,5a,5.,5.f (4.8) 

Y=[p, rf (4.9) 

U = [5ac,5rcf (4.10) 


The various variables and stability derivatives used here are explained in 
Appendix 1. \i/ and (}) are integrals of yaw rate and roll rate, respectively. 

The lateral (roll-yaw) autopilot is required to command stability axis roll rate 
and minimize induced sideslip. The stability axis co-ordinate system is defined as 
the transformation of the body axes to the stability axes using angle of attack (a) 
[7]. Thus roll rate in stability axis (ps) is given by 

Ps = pcos(a) + /■sin(a) . (4.11) 


4.3 Control Law Based on Decoupling and Pole-Zero Assignment 

Consider the multivariable feedback system shown in Fig.4.1, where P and 
C e C" ^ ” denote the plant and controller transfer function matrices, respectively. 
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r C(s) “ . P(s) 


Figure 4.1 : Multivariable Feedback Control System 

Let the plant of the system be factored as 

Pis) = Bis) = , (4.12) 

where, the pairs (A(s) , Bis)) and (Biis) , A 1 ( 5 ')) are the left and right coprime 
polynomial matrix factorizations of the plant, respectively. 

The sensitivity matrix, with the loop broken at the output of the plant, is given 

by, 

Sis) = iI + Pis)Cis))-^ (4.13) 

and then, the tracking error, Eis), is given by, 

Eis) = Sis) Ris) (4.14) 

where, Ris) is the reference signal. Let the reference signal be 

Ris) = [ri(5') , r2(5') ,.... , rn(j')] and let the zeros of the polynomial miis) be the 
poles of uis) in 7?e[5] > 0, for i = 1, ... , n. 

For the input-output decoupling, reference signal tracking and desired closed-loop 
pole assignment, the sensitivity matrix must be of the form. 

Sis) = diag ^ ii(5) , .... , snis) J 

_ j. >vi(j)/ni(i) Wnis)mnis) 

^i(^) ””” 5«(^) _ (4.15) 

where, giis), i = are Hurwitz polynomials with desired closed-loop 

poles and wiis), i = ,n, are undetermined polynomials which should be 

determined to satisfy the internal stability constraints. 

In the following. Lemma 1 and Lemma 2 due to [61] , are stated without 
proof. These provide pole-zero assignment control law which satisfies internal stability 
of the multivariable feedback system. . 


V 
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Lemma 1 

Suppose det (A(5')) and det (5i(j’)) have no common zeros in Re[s] > 0. Then 
S(s) is internally stable if and only if S(s) is analytical in Re[s] > 0 and for some 
appropriately dimensioned stable rational matrices X(5) and Y(s) 

S{s) = (4.16) 

I -Sis) = Biis)Xis) . (4.17) 

From equations (4.16) and (4.17), we have 

y(j) = 5(A A“ ^(5) (4.18) 

X(^) = (/-5(A) (4.19) 

Equivalently, S(s) is internally stable if and only if S{s), S{s) A~ ^( 5 ) and 
51 ^ (5) (/ - S(s)) are analytic in > 0. . 

REMARK : If S{s) satisfies the internal stability requirements, then the controller 
can be directly obtained from (4.13) as 

C(A = P~\s) sr\s) (I-S(s)) . . (4.20) 

Assume, 

aii(A ain(s) 

A-\s) = ■ ■ 

_a«l(A a„„(s) ^ (4.21) 

and 

*11 (A /'ln(A 

Bi~ \s) = • • . 

_ bnl(s) W . (4.22) 


From equations (4.15), (4.21) and (4.22), we obtain 
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au(^) JlM ain(5) 


S(J)A = 


Snis) anl(s) Snis) a,tn(s) 


(4.23) 


and, 


fcll(j) (l-i„(5)) bln(s) 

BT\s)(I-S(s)) = ■ ■ 

(l-si(s)) bnl(s) (1-Sn(s)) bnn(s) (4.24) 

Define the polynomials 

yi{s) = n/=i (s-pu) i = (4.25) 

where, Xi is the number of distinct poles pu of the ith row of A~ ^(s) in 
Re[ 5 ] > 0 and an is the greatest multiplicity of each pole pu which appears in any 

element of the ith row of A~ ^(s'). Similarily, we define the polynomials. 

Qiis) - nfL 1 (s - quf ^ , (4.26) 

where S/ is the number of distinct poles qu of the ith column Bj ^(j) in 
Re[5] > 0 and r\ii is the greatest multiplicity of each pole qu which appears in any 

element of the ith column of Bi ^(^). Then following Lemma is obtained : 


Lemma 2 

For the system shown is Fig.4.1, if yi(s) and ^i(s) are coprime for 
i = 1 , ... , n, then S{s) = diag [( 51 ( 5 ) , , 5'7j(^)] is internally stable if and only if the 
following conditions hold ; 

(i) Si{s) is analytic in Re [ 5 ] > 0, for 1 = 1, ... ,n; 

(ii) the numerator polynomial of si(s) is divisible by yi{s), for i = 1 , ... .n; 

(iii) the numerator polynomial of 1 - si{s) is divisible by Q.i{s) , 
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for I = 

remark 


If there exists any pair (yiis ) , Qi(s)) which is not coprime, then it is impossible 
for S(s) to achieve the internal stability. 


From the condition (ii) of Lemma 2, the numerator polynomial of Si(s) must 
contain yi(s), and from the equation (4.15), the numerator of si(s) must also contain 
miis). Thus the -numerator of si(s) must contain the least common multiple' Of yi(s) 
and mi{s) for i = I , .n ie. 


Sis) = diag 


Wl(^) OTl(j) 


Wn(s) mnis) 
8nis) 


his)ziis) lnis)Znis) 

s.W 


(4.27) 


where zi(s), i = 1, .... , n are the least common multiples of yiis) and mi{s) while 
/(•(s), i = 1, .... , n are undetermined polynomials. To satisfy the requirement of 
causality of the closed-loop system, the sensitivity matrix must be proper, i.e. 

deg igiis)) > deg {Uis)) + deg iziis)) . (4.28) 


From equation (4.27) we have 

Uis) = 1 - siis) 

^ 8i(s) - liis) Ziis) 

8iis) (4.29) 

and from condition (iii) of Lemma 2, the numerator of 1 - siis) must contain 
£2/(5). Thus we have 

hiis) = 8ii^) - Iiis)ziis) = Qiis)fiis), i = 1 .... ./i (4.30) 

where fiis), i = !,....« are undetermined polynomials. 

The following Lemma, given as Theorem 1 in [61], provides the condition on 
existence of Uis). 


Lemma 3 

The solution of Uis) in equation (4.30) exists if and only if nuis) is coprime 
with £2/(5), for i = 1 , ...n, respectively 

remark 



SYNTHESIS OF DECOUPLED ... ANNEALING 


73 


If li{s) exists and the number of undetermined parameters of li(s) is equal to 
deg(Qi(5)), then the solution of li{s) in equation (4.30) is unique for (Qi( 5 )), for 
i = 1, .... , n, respectively. 

By solving equation (4.30), we obtain li(s), i.e. S{s) is obtained. Then the 
controller can be derived using equation (4.20). 

Moreover, if li(s) exists and the number of undetermined parameters of li(s) is 
greater than deg(i2((5)), then no unique solution of li(s) exists for f = 1, ...,n, 
respectively. This leads to an over parametrized solution and the free parameters of 
li{s) can be determined according to some specific performance criterion. 

In the design here, li{s), fi{s) and gi{s) are obtained by solving nonlinear 
optimization problem formulated in Section 4.5 for the desired performance requirements. 


4.4 Performance Specifications 

In the following, the time domain and frequency domain performance specifications 
are expressed by performance functions ({)/ (/?), where R is the solution vector for 
the optimization problem. 

4.4.1 Time Domain Performance Specifications 

(i) Command Tracking 

Good tracking performance can be achieved by minimising the error between 

desired response and actual response. Let yfj{t) and yij(t) be the discretized desired 
and actual responses, respectively, for the ith channel for the discrete point j, 
j = 1, .... n?i. Then the performance function for command tracking is defined as : 

In order to get desired overshoot, another performance function (})2(7?) is defined 

for the error in desired overshoot (yf.max) and the actual overshoot (yi.max) and is 
given by 
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<t>2W — (yi , max ~ , max) • (4.32) 

The corresponding performance constraints are 

^i(^) = 0 (4.33) 

<^2iR) = 0 (4,34) 

(ii) Control Input 

The L 2 -nonn of the control input signal ui(t) for the ith channel, is minimized 
to save the control energy and to avoid saturation problems. 


The L 2 -norm of ui{t) is defined as [32] 

^ -, 1/2 


ll«:i|2 = 


r 

I u} (r) dt 




J 


(4.35) 


This L 2 -norm can be calculated using Parseval theorem [39] as given by the 
equation (2.17) of the Chapter 2. The corresponding performance function is defined 
as 


(}.3(i?) = (iiM,«li2-t/iol , (4.36) 

where C/jO is a constant which puts upper limit on L 2 -norm of uif) for the 
ith input-output channel of the closed loop system. The corresponding performance 
constraint is 

<i>3 {P) ^ 0. (4.37) 


(iii) Sideslip 

The sideslip is kept small to avoid large roll moments. The performance function 
for the sideslip is defined as 

(4.38) 

where hs is a constant which puts upper limit on the magnitude of permissible 
sideslip . tf is the time during which the sideslip is to be kept small. The 
corresponding performance constraint is given by 

<!>4 {R) < 0. 


(4.39) 
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4.4.2 Frequency Domain Performance Specifications 

(i) Robust Stability 

Robust stability specification is based on M-circle constraint [51]. 

The details of the M-circle constraint are explained in the Section 3.4.2 of the 
Chapter 3. The M-circle constraint, based on loop margin, imposes gain and phase 
margin restrictions in terms of M-circle radius. This constraint provides robust stability 
against plant parameter variations and time-varying nonlinearities. 

The performance specification, based on M-circle constraint for robust stability 
is defined as : 

10, ~) I 

where, Mmin is the minimum radius of the M-circle. 

The corresponding performance constraint is given by 
<^5{R) < 0 . 

(ii) Sensitivity and Complementary Sensitivity Shapings 
It has been shown in [48] and discussed in Section 3.3.2 of the Chapter 3 that 

many design specifications, time domain as well as frequency domain, can be achieved 
by proper shaping of the sensitivity and complementary sensitivity functions. 

Here, the sensitivity and complementary sensitivity functions at the output of 
the plant have been taken for the shaping. The performance requirements are the 
same as discussed in Section 3.4.2 of the Chapter 3. These are restated as follows; 

The spectmm of reference signals and disturbance signals are concentrated at 
low frequencies, while the spectmm' of the measurement noises extends over a much 
wider frequency range. Therefore, |S (/o))! is kept small at low frequencies 
(for 0 < 0) < coo , say) and \T (/co)| is kept small at high frequencies (for co > co^,, 
say). 

In a very low frequency interval 0 < co < coo (where coo is much lesser than 
the cut-off frequency), we require that : 

I Soi (/co)| < pi, for 0 < CO < coo 


(4.40) 


(4.41) 


and 


(4.42) 
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|7oi 0'ct>) 1 = 1, for 0 < CO < (00 . 


(4.43) 


Here, Soi (/(O) and Toi (/co) are the sensitivity and complementary sensitivity 
functions, respectively, at the plant output for the ith channel, pi (« 1) is a constant 
and puts upper bound on Soi (/co) in the frequency range 0 < co < coo- 

We define the performance functions ^eiR) and (j)? (R) corresponding to equations 
(4.42) and (4.43) respectively, as 




4)7 (R) = 


A Sup 


(0 e [0, cool 


{Toi (/co) -pi} 


(4.44) 

(4.45) 


Here, pi is a constant which puts upper bound on |roi(/co)| required in the 
frequency range (0, coo). 

The value of |ro: (/co)| at the bandwidth (co/,) is given by 

|ro: (j (Ob)\ = |ro/ (0)1/V2 = 0.707 ( 4 . 46 ) 

and the corresponding peformance function is defined as 

4)8 (R) = (\Toi (j cofc) 1 - 0.707/ (4.47) 

The loop transfer function. Lot, with the loop broken at the output of the plant, 
is P(s)C{s). The performance function corresponding to the value of loop transfer 
function at crossover frequency, coc, is given by 

4)9(i?) = /|in,(/coc)l-l/ (4.48) 


The permissible peak of the |5oj(/co)| has been specified in the perfomance 
function ^siR) of the robust stability given by equation (4.40). The permissible peak 
of [TOi (/co)| is specified by another performance function. 

(4.49) 

where pt is the bound on |roz (/co)| in the frequency interval CD e [0, co/,]. 

Therefore, to capture the desired shapes of |5oi (/cd)| and |roi (/cd)1 at various 
frequencies, we require. 
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4)6 (R) ^ 0 

(4.50) 

<^7 (R) ^ 0 

(4.51) 

bs (R) = 0 

(4.52) 

4)9 (R) = 0 

(4.56) 

4)10 (R) ^ 0 

(4.54) 


(iii) Disturbance and Measurement Noise Rejection 

As given by equation (A3. 6) in Appendix 3, the disturbance and measurement 
noise rejections can be obtained by keeping \Soi (/co)| small over lower frequencies 
and |roi 0'w)| small over higher frequencies, respectively. By proper sensitivity and 
complementary sensitivity shapings, both the disturbance and measurement noise 
rejections are achieved. However, explicit bounds on |Soi(/co)| and lroj( 7 ( 0 )| need to 
be put. 

Good disturbance rejection over the bandwidth of the feedback system can be 
achieved by keeping \Soi (/Ci3)| small over the bandwidth (0, coc). 


Hence, we define the following performance function for achieving disturbance 
rejection. 


A 

bii (R) = 


Sup 

CO e [0, =o] 


{\Si(j(p)\-bj((py , 


(4.55) 


where b/(CL)) is a continuous bound. The corresponding performance constraint 
is 

<^11 (R) ^ 0. (4.56) 


The choice of bf{(33) has been discussed in Section 3.4.2 of Chapter 3. The 
bound b/((0) is given by equations (3.43) and (3.44) of Chapter 3. For the sake of 
completeness of the specification, these equations are given below ; 

0 < 6/ = i/i « 1, if CD < COc 

bf = bfi > l, if CO > coc 

For the measurement noise rejection, the peformance specification is given by 

(4 5 ,) 

where c/ is a constant which puts the upper bound on \Toi (/a))| for frequencies 
0 ) > db- 


The corresponding performance constraint is given by 
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$12 (R) < 0 . 


(4.58) 


4.5 Formulation of the Design Problem 
4.5.1 Objective Function 

The objective function used, is the weighted sum objective function, same as 
was used in Section 3.4.3 of Chapter 3. It is given by the equation (3.70) of Chapter 
3. Here, the performance functions governing the inequalities are (j)3 (R), 4>4 (R), 
(j )5 (/?), ^6(R)^ (R)^ 010 (-K), 011 (I?) and 012 (.K), and the performance functions 

governing the equalities are given by 0i (R), 02 (R), 08 (R) and 09 (R) . 


4.5.2 Simulated Annealing (SA) Technique 

As discussed in the Section 3.4.4 of the Chapter 3, the SA technique guarantees 
global optimal or near global optimal solution. The same SA algorithm as was used 
in Chapter 3, is applied to solve the optimization problem, defined in Section 4.5.1. 


4.6 Design Methodology 

The closed-loop design goals for the roll-yaw MIMO system are to command 
stability axis roll rate (ps) and to minimize the side slip. The stability axis roU rate 
is given by equation (4.11). The control law discussed in the Section 4.3 provides 
input-output decoupling of the MIMO system. This facilitates independent tracking 
of body axis roll rate (p) and yaw rate (r). As given by equation (4.11), the command 
of the stability axis roll rate implies tracking of the body axis roll rate and yaw 
rate, i.e. 

Psc = Pc cos (a) + rc sin (a) (4.59) 

where psc, pc and rc are stability axis roll rate command, body axis roll rate 
command and yaw rate command respectively. In the present scheme of the autopilot 
design, p and r are used for feedback. Hence, the reference signal (r) here, is 
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and the measured output is ^ . Therefore, the roll-yaw missile dynamics form 
rc _ ^ _ 

a MIMO system with two inputs and two outputs. 

With input-output decoupling, the closed loop transfer function matrix is diagonal. 
The diagonal entries are the closed loop transfer functions for the individual channels 
and are given by equation (4.29). 

The., sensitivity function matrix, S{s), given by equation (4.27) is also diagonal. 

Therefore, the performance specifications, both time domain and frequency 
domain, can be achieved for each channel independently by proper selection of 
si{s) and ti{s). The pole-zero assignment is used to get proper si{s) and ti(s). As is 
evident from equations (4.27) and (4.29), the choice of poles and zeros for both 
si{s) and ti{s) is dependent on the choice of the polynomials U (s) and gi (s). The 
polynomial li(s), in turn, is related with the polynomial fi(s) by equation (4.30). 
Therefore, pole-zero assignment problem boils down to the selection of gi{s), li{s) 
and fi{s) to achieve the desired design specifications. 

The performance specifications are written as time domain and frequency domain 
specifications and the problem is posed as optimization problem. The objective function 
of the optimization problem is given by equation (3.70). of Chapter 3. 

This optimization problem, which is nonlinear, nonconvex and nondifferentiable 
in nature, is solved by SA technique to get gi(s), li{s) and fi(s)- The solution vector 
of the optimization problem consists of the unknown coefficients of polynomials 
§i(s), li(s) md fi(s). 


4.7 Design Algorithm 

Following are the steps in the design algorithm: 

Step 1 

Perform the factorization P{s) = A~ ^( 5 ) B(s) = 5i(.y)Al^(5) , as given in 
equation (4.12) and calculate A” ^(.y) and 51^(5'). Then determine yi{s) and Q.i(s) as 
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given by equations (4.25) and (4.26), respectively. Determine mi{s) and Zi{s) as in 
equations (4.15) and (4.27), respectively. 


Step 2 

Identify first order transfer function or second order transfer function with two 
complex conjugate poles and no zeros, to get desired unit step responses for body 

axis roll rate (p) and. yaw rate (r). Let these responses be yf (t) and yi(t), respectively. 

Step 3 

Specify the closed loop time domain and frequency domain performance 
specifications, including hard bounds, as given by equations (4.31)-(4.58). 

Step 4 

Using the performance functions as specified in Step 3, form the objective 
function with constraints, as given by equation (3.70). 

Step 5 

Set SA parameters and initial guess for the solution vector R. Choose initial 
values of non-negative scalar weight factors depending on the importance of design 
specifications. 

Step 6 

Run the SA algorithm to get optimal solution vector, R, consisting of unknov.-n 
coefficients of gi(s), li{s) and fi(s). 

Step 7 

Evaluate the performance with the solution as obtained in Step 5. If satisfactory', 
go to Step 9; otherwise go to Step 8. 

Step 8 

Modify the bounds on time domain and frequency domain performance 
specifications, and also modify the non-negative scalar weight factors depending on 
the accuracy desired. Reset SA parameters, if required. Go to Step 6. 

Step 9 

Obtain the controller, C( 5 '), as given by equation (4.20). 
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4.8 Design Examples 

Two design examples are worked out. The Example 1 relates with the BIT 
missile [7,62] for which both the longitudinal and lateral autopilots are designed in 
this thesis. Example 2 is for the HAVE DASH II [17], which is also a BTT missile, 
and is worked out to show the applicability of the design algorithm to other BTT 
missiles. 

4.8.1 Example 1 

The design example takes the roll-yaw missile dynamics described in Section 
4.2 of this Chapter. The flight conditions studied here are the same as considered 

for missile pitch dynamics, represented by a = 16*^, Mach 0.8 and an altitude of 4000 
ft. The nominal values of the dimensional aerodynamic stability derivatives for this 

flight condition, taken from [62], are reproduced below : 

Kp = - 0.0052 (1 A) 

Y5a = 0.1338(1 A) 

YSr = -0.1004(1 A) 

Lp = 556.8 (1 A^) 

L5a = 70.2(1 A) 

Lsr = -1079.0 (lA^) 

Vp = 5.679 (lA^) 

N5a = - 57.03 (lA^) 

Nsr = 5.7471 (lA^) 

The natural frequency and damping of the rudder and aileron, are taken the 
same as that of elevator, ie. (0=113 rad/s and ^ = 0.6. 

The transfer function matrix of the roll-yaw dynamics with the state space 
description of equations (4.6)-(4.10), is given by 

Pnis) 

P{s) = 

_P2l(s) P22(s) (4.60) 

where, 

P M = 8.9638e5(5|i4^ !e-l+2.0976el0(i + 5.3323e-l-2.0976el0 

' ’ 5ac s(s + 67 . st 4m-) (s + 67.8- 90.40 (s + 12. 1 69} (j - 12. 1 64) 
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and 


_ _2_ ^ -2.3993e7 (5 +1.7476e-2 + 2.66380 (j+ 1.7476e-2-2.663S0 
~ hrc~ s{s + 67.8 + 90.40 (^ + 67.8 - 90.40 (j + 12. 169) {s - 12. 164) 

_ j_ _ - 7.2822e5 {s + 12.462) {s - 12.47) 

~ hoc ~ i (i + 67.8 + 90.40 (j + 67.8 -90.40 (5 + 12.169) (s- 12.164) 


_ ji_ _ 7.3385e4 (s - 25.84) (s + 25.746) 

^22(^) - 5^^ - j (^ + 67.8) + 90.40 (5 + 67.8- 90.40 (i + 12.169) 12.164)' 


The plant transfer function matrix is unstable with one pole in RHP at 
12.164. 

Typical design specifications for a lateral autopilot can be stated as : 

Design a lateral autopilot to command stability axis roll rate with the steady 
state accuracy of 0.5%, time constant less than 0.2 sec, overshoot less than 20%. 
The sideslip should not be more than 5 degrees. The rudder and aileron should have 
deflections less than 60 degrees and fin rates should not exceed 400 deg/s. In each 
channel the autopilot must provide the minimum gain margins of (-6,+10)dB and 
the minimum phase margins of ± 40 degrees. The autopilot bandwidth is to be 
restricted by the high frequency unmodeled dynamics constraints, e.g. unmodeled 
flexible modes and actuator nonlinearities. 

The lateral autopilot is synthesised using the design algorithm given in Section 
4.6. 


Performing the factorizations of P(s), as in equation (4.2), we get 
ail (s) an{s) 

A-\s) = 


= 


a2l(j) a22is) 
bu(s) bnis) 


ms) 


b2lis) b22(s) 


and 


Ai(^) = 


1 

Ddis) 


dn(.s) 

421(5) 


412(5) 

422(5) 


(4.61) 

(4.62) 


(4.63) 


aii(5) = (5 + 6.3452) (5 + 4.2883) (5 + 3.279e-l + 1.8396) 

(5 + 3.2796e-l - 1.8396) (5 + 2.25) 
ai2(5) = 2.4439e2 (j + 6.0918) (5 + 3.9701) (5 - 2.0071) (5 - 2.1524) 
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021 W = 2.2054e-2 {s + 6.7987el) (j + 5.8929) {s + 3.8432) {s + 1.8052) 
a22{s) = (5 - 1 1. 1930) + 13.674) (5 + 5.9762) (^ + 4.0148) (^+ 1.9937) 

b\\{s) = 8.9838e5 (i + 5.8916) (i + 3.9556) (s + 2.0848) (5 + 51.565) 
bnis) = - 2.3993e7 (j + 9.9317e-l)(j + 6.2841) (j + 2.816) (5 + 4.5923) 
b 2 \{s) = - l.lilleS (s + 6.0829) (s + 4.027) (s + 1.9577) (5 - 1.5996e-l) 
b 22 (s) = 7.3385e4 (s - 3.7299e2) (s + 5.7005) (s + 3.4262) (s + 1.033) 
ifil(j) = (5 + 67.8 + 90.4i) (s + 67.8 - 90.4/) (s + 6.0741) (s + 4.022) (s + 1.9736) 

(j-7.3366e2 -3.0956/) ( 5 -7.3366e-2 + 3.09560 
di 2 (s) = 3.8175el (s + 67.8 + 90.4/) (5 + 67.8 - 90.40 (j + 5.7341) (j + 3.513) 

(s + 1.099) (s - 1.9576) 

= 2.0049e-2 (s + 1.1118c3) (s + 67.8 + 90.4/) ( + 67.8 - 90.40 + 6.0477) 

(i +4.0142) (i+ 1.9802) 

and 

d 22 (s) = (s + 67.8 + 90.4/) (s + 67.8 - 90.4/) (s - 10.173)(s + 5.8298) (5 + 3.6546) (s + 1.1314) 

where, 

Dd(s) = D(s) = (s + 1) (s + 2) (s + 3) (s + 4) (s + 5) (s + 6) (s + 7) 

Da(s) = sis + 67. S + 90.4/) (s + 67.8 - 90.4/) (s + 12. 169) (i + 12. 169) (^ + 12. 164) £)(j) 
Db(s) = 1.7406el3(^+1.8054e-2)r'(^) 

From the matrices A~ ^( 5 ) and 51 ^(s) we get the following values of polynomials 
Yl('y), Y2('S'), ^^l('S), and Q.2(s), using the equations (4.25) and (4.26): 


71 (-y) = Y2(^) = j(5- 12.164) (4.64) 

Qi(j) = Q2is) = 1 (4.65) 

Hence, we get 

miis) = m2{s) = s (4.66) 

and, 

21 ( 5 ) = Z2(s) = z(s) = 5(5-12.164). (4.67) 

The second order closed loop transfer function of provides the 


5*^ + IO 5 + 60 

desired time response. The closed loop transfer functions for both roll rate and yaw 
rate channels are taken of 5th order, same as that of the plant, satisfying the inequality 
(2.18) for the poles and zeros of the closed loop and open-loop systems. For the 
controller given by the equation (4.20) to be proper, the closed loop systems in both 
roll rate and yaw rate channels should have two zeros. 
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Here, we have = Q.2is) = 1 . Therefore, to satisfy the conditions (ii) and 

(iii) of Lemma 2, the requirements of proper controller and causality of the closed 
loop system as given in equations (4.20) and (4.28) respectively, we choose Zi (s). 


h (s), fl (s) and /2 (^) in the following forms: 


3 2 

ll(s) = s +a2S +aiJ + ao 

(4.68) 

/ 2 (j) = + a ' 2 + a'l s + a'o 

(4.69) 

fl(s) = b 2 S^ + bis + bo 

(4.70) 

and 


f 2 is) = b '2 + b'l s + b'o 

(4.71) 

Let gi (s) and g2 (s) be given by 


.gi (j) = / + ^4 + ^2 -y + ^0 

(4.72) 

g 2 (s) = + g'4S^ + g'zs^ + g'2S^ + g'lS + g'o 

(4.73) 


The gi(5), g 2 (s), li(s), his), fiis) and fzis) are related by the equality given in 
equation (4.30). The cofficients of all these polynomials are unknown. If all these 
coefficients are taken as design variables, in the optimization problem certain equality 
constraints among at, a'i, bi, b'i, gi and g'i will be involved to satisfy the equation 
(4.30). To avoid this problem, the design variables have been taken as au a'i, bi 
and b'i. Then, gi and g'i are calculated therefrom, putting the restriction that the 
poles of gi ( 5 ) and g2 (^) should lie in the LHP. Therefore, the solution vector, R, 
of the optimization problem is given by, 

R = [ai , a'i , bi ,b'if, V i = 1 , ... , 2 (4.74) 


The gain and phase margin constraints are imposed by M-circle constraint given 
by equation (4.39). Mmin = 0.7 guarantees gain margins of (-4.6, -f-10.4)dB and the 
phase margins of ±41 degrees, approximately. Initially, the design was aimed at 
achieving (Dc = 25 rad/s and CDfc= 45 rad/s. 

The following initial bounds were taken: 

Uio = 1.0 
U 20 =1.0 
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bs = 1.0 

bf\ = 0.4, for 0) < 0)c 
bf2 =1.4 

pi = 0.2 , for COO = 5 rad/s 

p2 = 1.0 

and ft 1.5 
cf = 0.7 . 

The following parametexs were set to run the Simulated Annealing algorithm, given in 
Appendix 4 : 

Ns = 20 
Nt = 50 

Ci - 2 , / = 1 , .... , n 

Nt = 4.0 

and, 

rr = 0.85 

The SA algorithm was run for different initial guesses. Initially, for the pole 
locations, one complex conjugate pole pair was taken as (s+5±5.9161i), which are 
the poles of the desired second order transfer function as in Step 2 of the design 
algorithm. The other poles were taken so that the above pole pair was in the dominant 
mode as compared to the remaining poles. Initial values of at, a'u bi and b'i were 
chosen to satisfy the equation (4.30). The bounds were set on au a'i, bi and b'i to 
have only positive real values. 

Initially, a set of scalar weight factors was selected and convergence of the 
SA algorithm was studied. The SA algorithm was run for various values of the 
temperature^parameter. As the starting temperature was increased, the solution improved. 
But, with increase of starting temperature speed of execution reduced. The starting 
temperature of 2.5el2 gave optimal solution within reasonable execution time. Only 
one run was required for the converged solution. • 

After each run, the performance of the design was evaluated with solution 
vector. And scalar weight factors and bounds were tuned, if required. The scalar 
weight factors were modified. 
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The scalar weight factors Xi, X 2 , X 3 , X4, ^5, ^5, Ay, As, A9, Aio, An and 
Ai 2 with the respective values of 1, 50, 100, 10, 100, 200, 200, 150, 30, 50, 42 
and 50 gave satisfactory performance for both the body roll rate and yaw rate 
channels. Finally, the solution converged with the bounds : 

Uio = 1.0 

U 20 = 0.5 


bs = 1.0 
bfi = 0.55 
bf2 =1.4 

pi = 0.4, for (00 = 10 rad/s 

pi = 1.0 

pi = 1.4 

and, 

cf = 0.7 . 

The values of the co£, and cOc were to be adjusted to 52 and 30 rad/s for both 
the roll rate and yaw rate channels . 


The following polynomials were obtained with the solution vector R: 
his) = s^+ 1.2314^2 s- + 2.7633e5 j + 1.4508e6 
his) = s^ + 8.3684el + 2..7897e5 5 + 1.0592e6 

Ms) = 1.5885e7j^+ I.1742e8 5 + 2.8716e8 
Ms) = 8.7777e5 /+ 1.2409e8 + 1.1926el0 1.9908ell 5 + 5.3544ell 

gl is)= / + 1.1098e2 5 ^^ + 2.7483c5 + 1.3975e7 + 9.9112el s + 2.8716e8 


and, 

g2 (i) = 5 ^ + 7.152el / + 2.7795e5 s^ + 1.2662e7 / + 8.049e7 s + 3.851e8 


(4.75) 

(4.76) 

(4.77) 

(4.78) 

(4.79) 

(4.80) 


The coiresponding controller, C{s), is given by the equation (4.20), as 
C(j) = p-J(j) [Sis)-I\ 

= [Si(^)AT'(5)r'[/-.r‘(5)] 


= Ai(^)PT'(j)[/-r*(^)] 


= 4i(j) Bi ^ (^) 


1 0 
0 1 


^l(.j) 0 

his)z\is) glis) 

0 his)z 2 is) 
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= Ai(5) by ' ( s) 


1 0 
0 1 


S(s) 0 
li(s)z(s) _ j(s )_ 

0 h(s)z(s) 


On substituting the values of Ai(5), Bi \s), gi (s), g2 (s), li(s), bis) and z(s) 
from equations (4.62), (4.63), (4.79), (4.80), (4.75), (4.76) and (4.67), respectively. 


we get 

Cu(-y) Ci2(.r) 

Dais) Dais) 

Cis) = 

Ciijs) Czzjs) 

Dais) Dais) (4.81) 

where. 


Cii(^) = 1.1657el2/+ 1.6658el4/+ 1.5283el6/ + 4.838el4i^- 1.0435el9 5^ 

-7.5124el9 5- 1.7901e20 

Ci2 is) = 3.598el4 / + 5.1042el6 / + 4.9116el8 / + 3.0393el9 + 1.5385e20 s^ 

- 1.741e20^-4.1495e20 

C 22 is) = 1.3442el3 / + 1.9208el5 / + 1.9129el7 / +2.1505el8 + 8.6317el9 

+ 4.9586e20i+1.9407e21 

C 21 is) = 1. 1568^13 / + 1.654el5 / + 1.577el7 / + 8.6185el7 - 2.2127el9 

-1.741e20 5-4.1495e20 

Del (i)= 1.7406el/+ 1.932el5^^+4.7838el8/-3.3168el9^^-3.0778e205--5.5458el85 

Dc 2 is) = 1.7406el3i® + 1.2452el5/ + 4.8381el/ - 4.0542el/ -2.25e20i^ - 4.0489el8j 

The unit step responses for p and r are shown is the Fig.4.2. Fig.4.3 shows 
the control deflections to unit step command and the sideslip is shown in Fig.4.4. 

Bode magnitude plots to indicate loop shapings near cOc for the loop transfer 
functions broken at the output channels of p and r are shown in Fig.4.5. Fig.4.6 
and Fig.4.7 show the gain characteristics of So (/co) and To (/co) at the output of p 
and r channels, respectively. The time domain and frequency domain performance 
specifications achieved with the designed autopilot are shown in Table 4.1. The 
maximum sideslip is 0.128 degrees. 


Dcfleclions (deg/) p.r (:rad/S) 
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Figure 4.2 ; Unit Step responses for Body Axis Roll Rate 
and Yaw Rate 
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Figure 4.3 : Control Deflections 
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Figure 4.4 : Sideslip 
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Figure 4.5 : Bode Magnitude Plots for Roll and Yaw Channels 
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Figure4.6 : Sensitivity and Complementary Sensitivity Functions for 
Roll Channel at the Output. 



Figure 4.7 : Sensitivity and Complementary Sensitivity Functions for 

Yaw Channel at the Output 
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Table 4.1 

Time Domain and Frequency Domain Performance for the Body Axis Roll 

Rate and Yaw Rate (Example 1) 


Channel 

Rise Time 

(s) 

Overshoor 

(%) 

Settling 
Time (s) 

Fin 

Deflection 

Fin Rate 
(deg/s) 

Gain 

Margin 

Phase 

Margin 

Gain 

Crossover 

Body 




(deg) 

(dB) 


Freq. (Hz) 

Axis Roll 

Rate 

0.04 

25.5 

0.4 

20.7 

469 

-11.3,+8.4 

± 74.3 

9.1 

Yaw Rate 

0.04 

31.8 

0.48 

1 1.42 

144 

I -10.6, +7.2 j 

± 73.3 

j 8.4 j 


Normally for the BTT missile, the maneuvering capability required in yaw is 
very small. The stability axis roll rate command tracking response can be obtained 
with combination of step responses of p and r using equation (4.11). If we take one 
to one combination of unit step responses for p and r, as shown in Fig.4.5. we get 
the unit step response for ps as shown in Fig. 4. 8. The design methodology facilitates 
any combination of p and r to get unit step response for ps. Thus, the present design 
methodology facilitates stability axis roll rate command tracking design for any 
capability in yaw. 



0.2 - 
0 


0 


' ' ! I I I 

0-5 I 1.5 2 2.5 3 

Time (sec) 

Figure 4.8 : Unit step Response for Stability Axis Roll Rate 
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Most of the design specifications have been reasonably achieved using present 
methodology of design except the overshoot. The RHP pole of the transfer function 
matrix of the roll yaw plant occurs in the polynomials 21 ( 5 ) and Z2(-S')- The zeros of 
the closed loop transfer functions for body axis roll rate and yaw rate channels are 
given by the roots of the polynomials fi(s) and fi(s), respectively. Due to RHP pole 
occuring in 21 ( 5 ) and Z2{s), the zeros of fi(s) and fzis) can not be assigned arbitrarily, 
as given by equation (4.30). As a result of which this RHP root forces the time 
responses in both p and r to be very fast and thus producing large overshoot. In 
both body axis roll rate and yaw rate unit step responses, the time constant is less 
than 0.1 sec and overshoot is more -than 25%. 

fiy varying the bounds on performance functions and scalar weight factors, 
different sets of design specifications can be achieved. 

4.8.2 Example 2 

In this Section, roll-yaw autopilot has been designed for the HAVE DASH n 
missile system . HAVE DASH n is a high performance BTT missile. The missile 
model has been taken from [17]. The linearized roll-yaw dynamics are described by 

p = Pj sin (a) - r cos (a) + fp- P + ( ^ + 76 a- 6a + kg • 6r 


Ps — Lp'P + Z/5 a • 5a + LS r • 5r (4.83) 

'r = N^-^+ ma-8a + N5r-5r (4.84) 

(? = Ps (4.85) 

8j = — CD (6a — 5ac) (4.86) 

6r =-CO(5r-5.c) (4.87) 


Here, ps is the stability axis roll rate ; 0o is the missile pitch attitude angle. 
The stability derivatives and other variables occuring in the above equations have 
already been explained with reference to roll-yaw missile model in Example 1. 

The missile is assumed to fly at 40,000 ft, Mach=2.75, a = 10^. At this flight 
condition the stability derivatives are as follows : 

7p = - 0.2980 (lA) 

75 a = 0.0025 (lA) 
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Ysr= 0.1021(1 A) 

Lp = -2075.6121 (1A-) 

Lsa = -2631.4756 (lA) 

Lsr = -486.4590 (lA^) 

Afp = 8.1743 (lA^) 

NSa = 4.2661 (lA^) 

N?,r = -117.2795 (lA^) 

A stability axis roll rate command of 90 deg/s has been taken between 2.0 to 
4.0 secs. The body yaw rate for coordination is 
rc = pjc'Sin (a) = 90 sin (10) = 15.62 

The ps and r are the measured variables and have been used for feedback. 
Therefore, the transfer function matrix of the roll-yaw dynamics , given by equation 
(4.60), is 

^12(^)" 

P22is)_ ( 4 . 88 ) 

-2.111e4 (j+ 1.1083e-l + 1.9139el i) (^ + 1.1083e-l - 1.9139el i) (j + 6.9224e-2) 
j (j+ 1.1514e-l + 1.9195el i) ( 5 + 1.1514e-l - 1.9195el i) (s + 6.7711e-2) ( 5 + 180) 

pn = 

Oac 

7.679e2 (j - 22.581) {s + 2.2814el) (s + 6.9218e-2) 

5 (i-+ 1.1514^-1 + 1.9195el z) (j+ 1.1514e-l - 1.9195€l z) (s + 6.7711e-2) (j+ 180) 



-8.7563e4 s{s + 3.6682e-l + 2.1311 e\ z) {s + 3.6682e-l - 2.2377el z) 
j (s + 1.1514e-l + 1.9195el z) (^ + 1.1514e-l - 1.9195^1 i) (j + 6.7711e-2) (j + 180) 


-4.1361e5s{s+ 1.4999e-l +2.1711 z) (j+ 1.4999e-l -2.1711 1) 
s (i + 1.1514e-l + 1.9195el z) (j + 1.1514e-l - 1.9195el i) (j + 6.7711e-2) ( 5 -+ 180) 


P{s) = 




where, 


pn = 
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The lateral autopilot is synthesised using the design algorithm given in Section 
4.6. Performing the factorizations of P(s) as in equation (4.2), we get 

ail (i) ai2(5) 

4 -^ 5 ) = 


and 


where, 


a2i(^) aiiis) 
bnis) bi 2 is) 




Db(.s) 


b 2 \{s) b 22 {s) 


A 1(5) 


1 


4ii(s) d\ 2 {s) 


Ddis) 


d 2 \is) d 22 {s) 

aii(5) = (5 + 180) {s + 4.9913) {s - 1.0658 + 1.8822el i) (j -1.0658 - 1.8822al i) 


(4.89) ’ 

(4.90) 

(4.91) 


(5 + 2.9624)(5 + 6.9853e-2) 

a\2(s) = -2.191 e-2 {s + 180) (s + 4.9915el) {s + 5.0049) {s + 2.996424) (s + 7.0434e-2) 

a 2 i(i) = 5.9443e2 ^ (i +180) (5 + 6.8887) (i + 4.3131) (i + 2.4318) 

a 22 (s) = sis + 180) is + 5.5126) is + 3.5144 + 1.294 i) is + 3.5144 - 1.294 i) (j + 2.8646) 

bnis) = -4.7367e5 i is + 5.4015) (5 + 3.6279) is + 1.8407) 

bnis) = -7.679e2(^ + 9.1862e2)(^ + 4.6042) (jr + 2.3435) (j + 3.7581e-l) 

b 2 lis) = - 8.7563e4 5 (5 + 5.2832) is + 3.4884) (i + 1.7445) 

b 22 is) = -2.11 le4 is + 4.7494) (j - 3.4097) is + 2.5327) is + 3.9942e-l) 

duis) = is + 7.195e-l + 2.267 I) is + 7.195e-l - 2.2673 t) is + 180) is - 4.33) is + 4.7335) 


dnis) = -3.3064el(j + 5.1521e-l+2.53120(j + 5.1521e-l-2.5312j)(5+180) 

(s + 4.5813) ( 5 + 2.2875) 


^^2l(^) = 1.0653 (j+ 180) (i- 1.1398el + 1.8453el z) (i + 1.1398el - 1.8453el 0 

( 5 +4.4779) ( 5 + 2.2405) 

and, 

d 22 is) = is + 180) is - 6.7915e-l + 1.9678el i) is - 6.7915e-l - 1.9678el i) is + 1.162el) 

(^ +4.4779) (s + 2.2405) 

where, 

Ddis) = Dis) = (j + l)(j + 2)(j + 3)(j + 4)(j + 5)(j + 6) 

Dais) = sis + 1.154e-l + 1.9195 I) is + 1.154e-l - 1.9195 i) is + 180) ^is+6.11e-3) Dis) 
Dbis) = 1.0066el0(^ + 2.9583e-l)D(j) 
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From the matrices A ^( 5 ) and Si ^(s), we get the following values of polynomials 


Yi(5-), yiis), and Qiis) using equations (4.25) and (4.26): 

Yi(i) = yiis) = i 
Qi(^) = Q.2(s) = 1 


and. 


Hence, we get 

mi (5) = m 2 {s) = s 

zi(5) = ziis) = z(s) = 5. 

The first order closed loop transfer function of 


20 


s + 20 


(4.92) 

(4.93) 

(4.94) 

(4.95) 

provides the desired time 


response. The closed loop transfer functions for both the roll rate and yaw rate 
channels are taken of 5th order, same as that of the plant, satisfying the inequality 
(2.18) for the poles and zeros of the closed loop and open loop systems. For the 
controller given by the equation (4.20), to be proper, the closed loop systems in both 
roll rate and yaw rate channels should have two zeros. 

Here, we have ^2i(5) = ^ 22 ( 5 ) = 1 . Therefore, to satisfy the conditions (ii) and 
(iii) of Lemma 2, the requirements of proper controller and causality of the closed 
loop system as given in equations (4.20) and (4.28), respectively, we choose h ( 5 ), 
h (s), fl (s) and fz ( 5 ') in the following forms; 

(4.96) 


and 


/l(j) = + az + a 2 S~ + a\ s ad 

his) = + a 2 + a '2 P" + a'is + a'o 

/l(-s) = b2S^ + b\s + bQ 


f 2 {s) = b '2 -^b'l s + b'o 


(4.97) 

(4.98) 

(4.99) 


The closed loop poles for both roll rate and yaw rate channels are chosen to 
be five simple poles . Then, ^1 ( 5 ) and gz (s) are of the following forms: 


^1 (-J) = (•y+Pl ){s + p2){s+p3){s+pi)is + p5 ) 

5 4 3 2 

= S +g4S +g2S +g2S +gis + g0 

glis) = (s+p'\)(s +p'2 )(s+p'3) {s + p'4)is+p'5) 


(4.100) 
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= ^ + ^' 35 ^ + g'2s^ + + s'o (4.101) 

The gi(^), g 2 {s), his), his), fiis) and fzis) are related by the equality given in 
equation (4.30). In this example one zero which appears in the denominator of 
P'^is) at -2.9583e-l has been taken in the polynomials /i ( 5 ) and fz (r) . Remaining 

one zero of fi is) "s and five simple poles of gi ( 5 ) 's have been taken as design 

variables. To avoid the introduction of additional constraints in the optimization 

problem due to equation (4.30), h ( 5 ) and h (r) have been calculated from the 
8 1 is), g2 is), fl ( 5 ) and fi (s) . The restriction has been imposed to have the roots 
of h is) and h is) in the L.H.R 

Therefore, the solution vector, R, of the optimization problem is given by, 

R = Ipi, p'i , Zi , z'i ]^, V / = 1 , ... , 5 (4.102) 

Here, zi and Zi are the zero locations in the polynomials /i ( 5 ) and fz is), 
respectively. Mmin of 0.7 guarantees desirable gain and phase margins, as discussed 
in Example 1. 

The initial values of bounds and scalar weight factors were taken same as in 
Example 1. Initially, the design was aimed at achieving (Oc = 25 rad/s and 0)^,= 45 
rad/s. 

The same set of SA parameters, as in Example 1, were taken. With initial 
temperature of 1.0e6, the optimal solution was obtained which gave satisfactory 
performance. 

The scalar weight factors Xi, X 2 , X 3 , A 4 , X 5 , Xe, X 7 , Xg, X 9 , Xio, Xn and 
ki 2 with the respective values of 1, 5, 500, 500, 70, 100, 100, 50, 30, 50, 42 and 
20 gave satisfactory performance for both the body axis roll rate and yaw rate 

channels. Finally, the solution converged with the bounds : 

U\Q = 1.0 deg 
U 20 = 2.0 deg 
bs = 2.0 
bfi = 0.25 
bfz =1.4 

p\ = 0.4, for coo = 10 rad/s 
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p2 = 1-0 
Pi = 1-4 

and, 

cf = 0.7 . 

The values of the Ob and coc were to be adjusted to 40 and 25 rad/s for both 
the roll rate and yaw rate channels . 

The following polynomials were obtained with the solution vector R; 


his) = / + 7.4009e2 + 2.0545e5 + 2.5354e7 i + 6.32 lle7 (4. 103) 

l2{s) = / + 2.6574e5 5^ + 2.3496el0 + 6.516el4 i + 1.8468eI4 (4. 104) 

/l(5) = 1.0759e4^^+ 1.114e95 + 3.2956e8 (4.105) 

fz(s) = 3.9826el3 s- + 1.3501el 5 + 3.9904el5 (4.106) 

gl (s)= / + 7.4009e2 / + 2.0545e5 + 2.5364e7 + 1.1772e9 j + 3.2956e8 (4. 107) 

and, 


g 2 (.s) = / + 2.6574e5/ + 2.3496el0r^ + 6.9143el4i^ + 1.3685el6j + 3.9904el5 (4.108) 
The corresponding controller, C{s), as in equation (4.81), is obtained with the 
following polynomials: 

Cii (j) = 5.0963e9/ + 5.286el4i^ + 9.5297el6/ + 5.9141el6 5^ + 4.5903el7.r^ 

+ 1.3308el7j-8.4553e-l 

Ci2 (s) = 3.0582el6 / + 1.5881el9 / + 1.8582^21 / - 7.0598e21 / - 9.6327e23 

-3.5078e23 5 - 1.9668e22 

C 21 ( 5 ) = -9.4212e8 - 9.7718^13 / - 1.766el6 / - 6.7039el6 - 8.8125el8 

-2.6016«18i + 8.7086 

C 22 (i) = 8.4073el7 / + 4.3658e20 / + 5.182e22 / + 1.8997e23 + 1.8838e25 

+ 6.8574e24j + 3.844e23 

Del (s) = 1.0066el0/H-7.4527el2.s^ + 2.07el5s'^ + 2.5582el7? + 7.1178el7j^+ 1.88el7i 
Del is) = 1.0066el05^ + 2.615el5s^ + 2.3652e20/ + 6.559e24j^ + 3.7994e245^ + 5.5_e23j 

The stability' axis roll rate response due to roll rate command of 90 deg/s, 
between 2.0 and 4.0 secs, is shown in Fig.4.9. Fig.4.10 shows the yaw rate response. 
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Figure 4.1 1 ; Aileron and Rudder Deflections 
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Figure 4.12 ; Step Response for Stability Axis Roll Rate 
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Figure 4.16(a) : Magnitude Plots for Stability Axis Roll Rate 

and Yaw Rate Channels 
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Figure 4.16(b) : Phase plots for Stability Axis Roll Rate Channel 

and Yaw Rate Channel 
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Figure 4.18 ; Sensitivity and Complementary Sensitivity Functions for the 
Yaw i^.ate Channel at the Output of the Plant 

stability axis roll rate due to 90 deg/s is shown in Fig.4.12. The sideslip 
response due to the initial condition is shown in the Fig.4.13. The Fig.4.14 
shows the sideslip due to stability axis roll rate command of 90 deg/s. Fig.4.15 
shows the sideslip due to roll rate command of 90 deg/s between 2.0 and 4.0 
secs. The maximum sideslip is 1.9 deg. 

Figs.4. 16(a) and 4.16(b) show the Bode magnitude and phase plots for the 
stability axis roll rate and yaw rate channels with the loop broken at the output of 
the missile plant. Figs.4. 17 and 4.18 show the gain characteristics of the sensitivity 
and complementary sensitivity functions for the roll rate and yaw rate channels, 
respectively, at the output of the plant. 

The time domain and frequency domain performance specifications achieved 
with the designed autopilot are shown in Table 4.2. The maximum sideslip is 1.9 
degrees. 

4.8.3 Results and Discussion 

Time domain and frequency domain performances for both the designs in 
Examples 1 and 2 are shown in the Tables 4.1 and 4.2, respectively. In Examplel, 
the design specifications have been reasonably achieved except overshoot, which is 
on higher side. High overshoot results due to limitation imposed by the RHP pole 
at 12.164 on time response and bandwidth. As pointed out by Middleton [76], the 
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Table 4.2 

Time Domain and Frequency Domain Performance for the Body Axis Roll 

Rate and Yaw Rate (Example 2) 


Channel 

Rise Time 

(s) 

Overshoot 

(%) 

Settling 
Time (s) 

Fin 

Deflection 

(deg) 

Fin Rate 
(deg/s) 

Gain 

Margin 

(dB) 

Phase 

Margin 

(deg) 

Gain 

Crossover 
Freq. (Hz) 

Stability 
Axis Roll 

Rate 

0.08 

0.0 

0.16 

4.0 

20.2 

- 6,400 

± 93.442 

3.3 

Yaw Rate 

0.08 

0.0 

0.16 

2.5 

10.6 

-12.5, +14.3 

± 71.884 

7.1 


Table 43 

Tradeoff Among Rise Time, Overshoot and Bandwidth 


tn sec 

overshoot {%) 

B.W. (rad/s) 

0.02 

14.37 

115.0 

0.03 

23.43 

76.67 

0.04 

34.0 

57.5 

0.06 

60.43 

38.33 
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overshoot, bandwidth, rise time and RHP pole location are restricted by the 
following two relations: 

overshoot > [{p- tr- 1) + 1 ] 

p- tr 

and, 

trm = 2.3. 

The second one is the empirical relation. The variation of the rise time, percentage 
overshoot and bandwidth are shown in the Table 4.3. At p=12.164, to have the rise 
time of 0.04 sec, the overshoot and bandwidth will be of the order of 34% and 
57.5 rad/s , respectively. If the rise time is reduced, the lesser overshoot will be 
obtained but at the same time, the bandwidth will increase and vice versa. To reduce 
the overshoot at the cost of increased bandwidth is disadvantageous as it may excite 
high frequency modes of the missile dynamics and hence the robustness of the closed 
loop system against the neglected and mismodeled dynamics will reduce. On the 
other hand, since other time response specifications are good, the higher overshoot 
will not have significant effect. 

As can be seen from the Table 4.2, the achieved performances in time and 
frequency domain in Example 2 are excellent. In this case, the overshoot does not 
occur. This is due to the fact that the design was attempted with real simple poles 
in LHP, and it worked out giving excellent results. The tracking response for stability 
axis- roll rate is fast with the rise time of 0.08 sec and settling time of 0.16 sec.. 
Yaw rate response is also very fast which is required for turn coordination. Both the 
roll rate and yaw rate tracking responses have been independently achieved. The 
maximum sideslip is 1.9 degrees. With different sets of bounds and scalar weight 
factors, the wide variety of specifications can be achieved. 


4.9 Conclusions 

The design methodology developed in this Chapter, presents a systematic 
procedure of the autopilot design for the coupled roll-yaw dynamics for the stability 
axis roll rate command tracking, while keeping the sideslip small. The input-output 
decoupling has been achieved which facilitates stability axis roll rate command tracking 
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design with independent design of body axis roll rate and yaw rate channels In 
other words, the stability axis roll rate command design can be achieved with varied 
capabilities in yaw maneuver. The design tradeoff has been achieved by vaivinv the 
bounds on the performance functions and scalar weight factors during the' ttesign 
process. Therefore, fte wide variety of design specifications can be achieved by 
suitably choosing the performance bounds and scalar weight factors. The SA algorithm 

proved to be very effective for optimization. Only one run of the algorittei was 
required to get the converged solutions. 


Chapter 5 

Robustness Analysis 

5.1 Introduction 

Primarily singular value and structured singular value analysis techniques 
[62,64-66] are the recent tools to evaluate the robustness of the feedback control 
systems. 

Doyle [65] and Lehtomaki [66] developed singular value based robustness tests 
modeling uncertainties as a full single block matrix structure. If uncertainties are 
really unstructured and in this form, then these tests are not conservative. If the 
uncertainties occur in the system in a structured way, then the robustness tests produce 
conservative estimates of stability robustness. 

In this Chapter the singular value based analysis techniques have been used to 
estimate the robustness of the autopilot designs in the Chapters 2, 3 and 4. The 
uncertainties modeled for the analysis include neglected actuator dynamics, mismodeled 
actuator dynamics and neglected sensor dynamics. These uncertainties are modeled 
as unstructured uncertainties. 

S 

5.2 Robustness Theory for Singular Value Based Analysis 

The singular value based analysis for the unstructured uncertainty is based on 
the robustness theorems which are derived from an application of the multivariable 
Nyquist theorem [68]. The robustness problem using the multivariable Nyquist technique 
is to determine to what extent the parameters in the loop transfer matrix (LTM) can 
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vary without compromising the stability of the closed loop system. It is shown in 
[68] that 

det [/ + L(s)] = (J)d (s)/ <1)01 {s) , (5 

where (!)cl (s) and ({)ol (s) are open loop and closed loop characteristic polynomials, 
respectively. 

The fundamental theorem of the robustness [62,66] is given below : 

Theorem 1 [62,66] 

The polynomial (^'d (^) has no closed right half plane zeros, and the perturbed 
feedback system is stable if the following holds : 

(a) 4 ) 0 / (s) and ^'ol ( 5 ) have the same number of closed right half plane zeros; 
(])ci (5) has no closed right half plane zeros. 

(b) det [/ + L(s , e)] = 0, V ( 5 , 6 ) in Dr x [0, 1] and V R sufficiently large, 

where ^'ol (s) and (s) are the polynomials ^ol (s) and ^d ( 5 ) with uncertainties. 
L{s , e) is the loop transfer function matrix with real coefficients which are continuous 
in 8 for all 8 such that 0 < e < 1. is Nyquist contour which encloses all closed 
right half plane zeros of <^ol (s). L(s , 8) which satisfies 

Lis, 0) = Lis) (5.2) 

and, 

Lis, 1) = L'is) (5.3). 

This theorem is used to develop the test for different types of model 
characterizations. The most common model error characterizations are additive errors 
and multiplicative errors. 

Let Eis\ denote the modeling error under consideration. The additive model 
error is given by 

Eais) = L'is) -Lis) (5.4) 

and, the multiplicative model error is given by 

Emis) = [L'(j) - L(j)] L- * is) (5.5) 

The perturbed LTM for additive errors is given by 
Lis, e) = L(s) + e Ea{s) 

and, the multiplicative error is given by 


(5.6) 
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L(s, e ) = [/ + e SnCj )] L(s) 

The following are the stability robustness theorems for robustness analysis due 
to Lehtomaki [62,66] 

Theorem 2 (for additive uncertainty model) 

The polynomial ^'d (s) has no closed right half plane zeros and the perturbed 
feedback system is stable if the following holds: 

(i) <t>c/ (•!?) has no closed right half plane zeros. 

(ii) a [/ + L(r)] > o [£«(.?)], V r e Dr and V R sufficiently large, with Eais) 
given by equation (5.4). a( . ) and g( . ) are smallest and largest singular values, 
respectively. 

Theorem 3 (for multiplicative uncertainty model) 

The polynomial (j)'c/ (s) has no closed right half plane zeros, and the perturbed 
feedback system is stable if the following holds: 

(i) ^dis) has no closed right half plane zeros. 

(ii) a [J + L" ^(.y)] > ^[^^(y)], V y e and V R sufficiently large, with 
Em{s) given by equation (5.5). 

The proofs of the theorems 1 and 2 are available in [62,66]. These theorems 
provide sufficient tests of stability. As long as the singular value frequency responses 
do not overlap, the stability is guaranteed. 

S 

5.3 Autopilot Robustness Analysis 

In this Section, the singular value based robustness theory is applied to the 
uncertainties due to neglected and mismodeled actuator dynamics, and neglected sensor 
dynamics. These uncertainties are modeled as unstructured multiplicative uncertainties. 
The autopilot sensitivity is analysed for the designs in Chapters 2, 3 and 4. 
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Figure 5.1 : Feedback Control System 

5.3.1 Actuator Uncertainty Analysis 

Here, the autopilot sensitivity is analysed for the neglected actuator dynamics 
and mismodeled actuator dynamics in respect of the designs in Chapters 2, 3 and 
4. 

5.3.1.1 Neglected Actuator Dynamics 

Fig.5.1 shows the autopilot transfer function C(^) and the uncertain missile 
dynamics G{s). For the multiplicative plant uncertainty at the plant input, we have 
G(i) = Go{s) [l+Em{s)\ , (5.8) 

where Gq{s) is the open loop transfer function matrix of the missile and Em{s) is 
the multiplicative uncertainty due to neglected actuator dynamics at the plant 
input. The loop transfer function at the plant input, for the true missile model is 
L\s) = C{s)G{s)[l + Em{s)] , ( 5 . 9 ) 

and for the nominal design model is 

L{s) = C{s) G{s). (5.10) 

Our aim is to study the robustness of the pitch and roll-yaw designs to neglected 
actuator dynamics. For this, the Theorem 3 of the Section 5.2 is applied. That is, 
the following relationship should hold for the stability robustness : 

a[/ + L-^5)] > a[£;„(5)], (5.11) 

with L{s) given by the equation (5.10). This singular value based test is only 
sufficient condition for stability. As long as the singular value frequency 

responses of a [7 + L” ^ ( 5 )] and a [Em ( 5 )] do not overlap, the stability is 
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guaranteed. The robustness test is applied here to find out the lowest actuator 
natural frequency of the actuator dynamics that guarantees stability. 

Now the robustness test is applied to the longitudinal autopilot designs in 
Chapters 2 and 3, and the lateral autopilot design in Chapter 4. 


5.3.1.1.1 Longitudinal Autopilot Design (Chapter 2) 

The actuator in pitch dynamics is the elevator which is modeled as the second 
order transfer function. The mutiplicative uncertainty created by neglecting the actuator 
dynamics is given by, 

r f \ - _ 1 _ - i (.y + 2 ^ (Oa) 

OM-S) / + 2^0)aJ + (DS (5.12) 

where ^ = 0.6 and (Da is taken as the unknown actuator natural frequency. Sg 
and 5ec are elevator fin deflection and commanded elevator fin deflection, 
respectively. Here, C(s)' = [Ci(s') C2(5) C2(5')] , where Ci (^) and C 2 (s) are given 
by equations (2.26) and, (2.33), and G(s) is given by equation (2.20), 
respectively, of the Chapter 2. The application of the equation (5.11) is shown in 
Fig.5.2. Three natural frequencies of actuator dynamics at (t)a=2, 3.2 and 5 rad/s 
are taken for the plots of a [E,n(s)]. The intersection of frequency response 

curves of o[I+L~^ ( 5 )] and a [Em ( 5 )] violates the condition given by equation 
(5.11). It is evident from the Fig.5.2 that for cOa = 3.2 rad/s, a[£'/n(5)] curve 

is just touching a[/ + L~^(5)] curve and for 0)a < 3.2 rad/s the robustness 
condition is violated. This shows that the stability of the missile pitch control 
system is not guaranteed for the actuator natural frequencies below 3.2 rad/s. 

N 

5.3.1. 1.2 Longitudinal Autopilot Design (Chapter 3) 

The application of the robustness test for the neglected actuator dynamics, given 
by equation (5.11) is shown in Fig.5.3. G{s) and C{s) used here for the robusmess 
analysis are the same as shown in Fig.3.2 of Chapter 3. a [Em (^)] is plotted for the 
four natural frequencies of actuator dynamics at cOa = 2, 3.2, 5.0 and 10 rad/s. For 
(Da=10 rad/s, the stability of the missile control system is not guaranteed. 
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Figure 5.2 : Pitch Robustness to Neglected Actuator Dynamics 

in Chapter 2 
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5.3.1.1-3 Lateral Autopilot Design (Chapter 4) 

The actuators in the lateral missile dynamics are aileron and rudder. The dynamics 
of the aileron and rudder is taken same as the elevator dynamics. Therefore, here 
the Emis) is a 2x2 matrix with diagonal elements given by equation (5.12). Now, 
the robustness of the lateral autopilot designs in Examples 1 and 2 of Sections 4.8.1 
and 4.8.2 is evaluated. 

(a) Example 1 

The application of the equation (5.11) for the neglected actuator dynamics is 
shown in Fig.5.4(a). a [^/^(s)] is plotted for three natural frequencies of actuator 
dynamics at (Ma = 20, 60 and 120 rad/s. It is clear from the Fig.5.4(a) that the 
missile control system in roll-yaw does not guarantee stability for the actuator natural 
frequencies below 20 rad/s. 

(b) Example 2 

Fig.5.4(b) shows the application of the equation (5.11) for the neglected actuator 
dynamics. It is clear from the Fig.5.4(b) that the missile control system in roll-yaw 
does not guarantee stability for the actuator natural frequencies below 60 rad/s. The 
missile model taken in this example has first order actuator with the natural frequency 
of 180 rad/s. Therefore, the robustness against the neglected actuator dynamics is 
good. 

5.3.1. 2 Mismodeled Actuator Dynamics 

5.3. 1.2.1 Longitudinal Autopilot Design (Chapter 2) 

The natural frequency of the actuator dynamics is mismodeled by the relation 
(Utrue = e (Oa [62]. Here, ooa is the nominal actuator frequency and is assumed to 
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be 67.8 rad/s. The multiplicative uncertainty model for the mismodeled dynamics 
is given by, 

r- 1(1 -e‘)i + 2^coae (1-e)] 

CmW - 2 - 2 2' 

S +2c,(i)aes + (ilat (5-13) 

The Fig.5.5 shows the application of the robustness test given by equation 
(5.11). a [i?/,! (5)] is plotted for four actuator frequencies at (£itrue= 3.14, 12.56, 31.4 
and 62.8 rad/s. The stability robustness condition given by equation (5.11) is violated 
for (x)true 3.14 rad/s. 


5.3. 1.2.2 Longitudinal Autopilot Design (Chapter 3) 

Fig.5.6 depicts the application of the condition given by equation (5.11). 
a[£'/«(^)] is plotted for the four actuator frequencies at untrue = 6.28, 12.56, 31.4 

and 62.8 rad/s. The missile control system in the pitch does not guarantee stability 
against mismodeled actuator dynamics for the actuator natural frequencies below 6.28 
rad/s. 

5.3.1.2.3 Lateral Autopilot Design (Chapter 4) 

(a) Example 1 

Here, the Em(.s) is a 2x2 matrix with the diagonal entries given by equation 
(5.13). Fig. 5. 7(a) displays the application of the stability robustness condition given 
by equation (5.11). o [E,n{s)] is plotted for four frequencies at (Utme = 10, 20, 60 
and 120 rad/s. The stability robustness condition is violated for (Utrue < 14 rad/s . 


(b) Example 2 

Fig.5.7(b) shows the application of the stability robustness condition given by 
equation (5.11). a [Em (5)] is plotted for the four frequencies at (Dtme = 10, 20, 60 
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This is good robustness against the mismodeled actuator, in this example, as the 
actuator model talcen in the missile dynamics is of the first order with the 
natural frequency of 180 rad/s. 


5.3.2 Sensor Uncertainty Analysis 

The uncertainty due to neglected sensor dynamics is modeled as multiplicative 
uncertainty at the output of the plant. 


5.3.2. 1 Longitudinal Autopilot Design (Chapter 2) 

The sensor dynamics has been neglected in the autopilot design. The accelerometer 
and the rate gyro are used for measuring the normal acceleration and pitch rate 
respectively. For the robustness analysis purpose, we model both the accelerometer 
and rate gyro dynamics using second order transfer function given by 
5s _ _ col 

Ssr .5-^+ 2 ^ 0 )^^ + col ’ (5.14) 

where 5 ^- and are sensor (accelerometer, rate gyro) deflection and commanded 
deflection respectively. ©5 and ^ are the sensor natural frequency and damping ratio 
respectively. ^ = 0.5 and ©5 is taken as unknown sensor frequency. The sensitivity 
of the autopilot design is analysed by finding out the minimum natural frequency of 
the sensor dynamics violating robustness condition given by equation (5.11). In 
equation (5.11), the L(s) becomes the loop transfer function with the loop broken at 
the output of the plant. 

Here, the error due to the neglected sensor dynamics, both for accelerometer 
and rate gyro, is given by 

-s{s + 2zs0^s) 

Einsx^f ■“ o 0 

i^ + 2^(Os5 + co| (5.15) 

The Em{s} in equation (5.11) becomes a 2x2 matrix with the diagonal entries 
given by equation (5.15). 
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Fig. 5. 8 shows the application of the equation (5.11) for the robustness analvsis 
due to the neglected sensor dynamics, o [£'/?j(5)] is plotted for three sensor natural 
frequencies at = 40, 80 and 200 rad/s. 

As can be seen from the plots, the plots of a [£'m(^)] for tUs = 40 and 80 

rad/s intersect with a[I + L~^ (^)] plot where as a [E'mC^)] plot for (0^= 200 rad/s 
does not intersect. It has been found that below =107.5 rad/s, the plots of 

o[I+L~^ (5)] and a [F,«(s)] always intersect with each other and hence the missile 
control system in pitch does not guarantee stability against neglected sensor dynamics 

below sensor natural frequencies of 107.5 rad/s. The ct [7 + L" ^ (5)] curve does not 

look smooth . It is due to the fact that in the computation of a [7 + L” ^ (5)], the 
poor condition number estimate was indicated. 

5.3.2.2 Longitudinal Autopilot Design (Chapter 3) 

Fig. 5.9 depicts application of the equation (5.11) for the robustness analysis 
against neglected sensor dynamics, a [^^(s)] is plotted for three sensor natural 
frequencies at ©j = 40, 80 and 200 rad/s. The plot of a [-E'/n(‘S')] for ©^ =40 and SO 

rad/s intersect with the o[I + L~^ (5)] plot whereas q [£'/n(5)] plot for ©^ = 200 rzdls 
does not intersect. It has been found that intersection always takes place between 

o [^^(s)] and c [7 + L~ ^ (s)] below ©5 = 102.6 rad/s. Hence, the robustness condition 
given by equation (5.11) is violated for ©^ = 102.6 rad/s. Again, here also, in ±e 

computation of the a [7 + L~ ^ (s)], the poor condition number estimate was indicated. 

5.3.23 Lateral Autopilot Design (Chapter 4) 

The roll-yaw autopilot design in Chapter 4 uses body axis roll rate and yaw 
rate for the feedback. The sensors employed are rate gyros for both the roll and 
yaw channels. The dynamics of the roll rate gyro and yaw rate gyro are taken same 
as given in equation (5.14). Now, the robustness against the neglected sensor dynamics 
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Figure5.8 : Pitch Robustness to Neglected Sensor Dynamics in Chapter 2. 
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Figure 5.9 ; Pitch Robustness to Neglected Sensor Dynamics in Chapter 3 
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Figure 5.11: Roll- Yaw Robustness to Neglected Sensor Dynamics 
for Example 2 
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is evaluated for the Examples 1 and 2 of the Sections 4.8.1 and 4.8.2 


(a) Example 1 

The application of the stability robustness condition given by equation (5.11) 
is shown in Fig.5.10. a is plotted for coj = 40, 80 and 200 rad/s. 

As is seen from the Fig.5.10, the a[£'/n( 5 ')] and a[7 + L~^( 5 )] plots intersect 
for (Oj=40 and 80 rad/s while they do not overlap for (0^=200 rad/s. It has been 
found that intersection always occurs for coj = 95 rad/s. Flence the lateral missile 
control system does not guarantee stability for the sensor frequencies below 95 rad/s. 


(b) Example 2 

The plots for a [/ + L”^ (s')] and a[Em(s)] are shown in Fig.5.11. a[£) 7 z(^)] is 
plotted for CDj = 40 and 80 rad/s. It is seen from the Fig.5.11 that ^[£^(■5')] curve 

at 0)5 = 80 rad/s intersects the g[/+L~^ (s)] curve whereas a [£’m('S’)] curve at 0)5=40 
rad/s just touches. Therefore, the stability robustness condition is violated for the 
sensor frequencies below 40 rad/s. 

5.3.3 Results and Discussions 

The results for the singular value based robustness analysis are given in Table 
5.1. It can be seen from the Table 5.1 that all the autopilot designs provide reasonably 
good stability robustness against neglected actuator dynamics and mismodeled actuator 
dynamics. The frequencies below which the robust stability is not guaranteed against 
the neglected actuator dynamics and mismodeled actuator dynamics for closed loop 
system in roll-yaw dynamics, for the Example 2 in the Section 4.8.2 of Chapter 4, 
are 60.0 and 55.2, respectively. Though frequencies look higher in comparision with 
the corresponding frequencies in other designs but the natural frequency of the actuator 
used here is 180 rad/s. Thus, these are also the good results for robust stability 
against neglected and mismodeled actuator dynamics. 
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Table 5.1 

Robustness Analysis Results 



Minimum Actuator 
Natural Frequency for 
Stability Against 
Neglected Actuator 
Dynamics (rad/s) 

Minimum Actuator 
Natural Frequency for 
Stability Against 
Mismodeled Actuator 
Dynamics (rad/s) 

Minimum Sensor Natural 
Frequency for Stability 
Against Neglected Sensor 
Dynamics (rad/s) 

Longitudinal Autopilot 
Design (Chapter 2) 

3.2 

3.14 

107.5 

Longitudinal Autopilot 
Design (Chapter 3) 

10.0 1 

6.28 

102.6 

Lateral Autopilot Design 
(Chapter 4. Example 1) 

20.0 : 

14.0 

1 

94.8 

I Lateral Autopilot Design 
[ (Chapter 4, Example 2) 

60.0 

55.2 

40.0 


The results against neglected sensor dynamics are relatively poor. For the 

longitudinal autopilot designs in Chapters 2 and 3, for keeping the closed loop system 

stable in pitch, the minimum sensor natural frequencies should be more than 107.5 

and 102.6 rad/s, respectively. In practice, the sensors are chosen with the natural 

frequencies higher than these values. For example, rate gyros selected for high 

performance missiles have natural frequencies of the order of 200 rad/s. But, still 

the margin available here is lesser than that for actuator natural frequencies. For the 

roll-yaw autopilot design, the minimum frequency for the sensor dynamics to guarantee 

\ 

stability are 95 rad/s and 40 rad/s for the Examples 1 and 2, respectively, which are 
reasonable. 

5.4 Conclusions 

The singular value based robustness analysis has been carried out to analyse 
the autopilot sensitivity to neglected and mismodeled actuator dynamics, and neglected 
sensor dynamics. These results, given in Table 5.1, show that all the autopilot 
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designs provide reasonably good stability robustness against neglected and mismodeled 
actuator dynamics. The results against neglected sensor dynamics are relatively poor 
for longitudinal autopilot designs in Chapters 2 and 3. 

The reason for the poor results for the stability robustness against the neglected 
sensor dynamics for the longitudinal autopilot designs can be attributed to the fact 
that frequency loop shaping in the longitudinal autopilot designs has been earned 
out at the input of the plant. For improving the stability robustness against neglected 
sensor dynamics, the frequency loop shaping should be carried out at the output of 
the missile plant also. 



Chapter 6 
Conclusions 


6.1 General 

The bank-to-tum missiles are high performance missiles employed for air to air 
and air to ground applications. The BTT configuration taken for autopilot design in 
this thesis work is meant for air to ground applications. These missiles have to 
perform under changing and uncertain environments. This poses a challenging task 
for the autopilot design. The objective of the work reported in this thesis has been 
the parametric synthesis of the autopilot for BTT missile. 

In this Chapter, the results obtained in the work reported in this thesis are 
reviewed and the contributions made are summarised. Suggestions for further scope 
of the work in this area are also indicated. 

6.2 Review of The Work Done 

In Chapter 2, the parametrized longitudinal autopilot has been designed using 
pole-zero assignment. The systematic design methodology and algorithms have been 
developed. The pitch dynamics form a SIMO system with one input and two measured 
outputs, namely normal acceleration and pitch rate. In the design methodology 
developed here, one loop has been designed at a time. The inner pitch rate loop has 
been designed for stabilization, damping and minimum control input whereas the 
outer acceleration loop has been designed for acceleration command tracking with 
specified gain and phase margins at the input of the missile plant. The control law 
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is based on the Youla parametrization of stabilizing controllers which facilitates 
incorporation of poles and zeros through the proper choice of the polynoniials. The 
poles and zeros are selected solving a nonlinear optimization problem. The objective 
function includes the terms corresponding to command tracking error, overshoot, 
undershoot, gain and phase margins. The objective function is non-differentiable and 
has been solved using the direct search methods of Simplex search, Powell’s conjugate 
direction and the Pattern search. The methodology developed, provides systematic 
synthesis procedure for the autopilot design in which the subjective judgement is 
minimized but the drawbacks are : (i) the order of the controller is high, and (ii) 
the conventional search methods used here, do not guarantee global optimum solution 
but only local optimum. 

In Chapter 3, the longitudinal autopilot has been designed using parametrization 
of stabilizing controllers in two degrees of freedom configuration. The possible classes 
of stabilizing controllers and the desired characteristics are provided by the choice 
of two free parameters. The tracking performance depends on one parameter whereas 
the feedback properties depend on the other parameter. The parameter appearing in 
the tracking performance has been selected using closed loop model matching concept 
for the desired time domain performance. The parameter appearing in the feedback 
properties has been selected for the desired frequency loop shaping. The frequency 
loop shaping has been done at the input of the plant. The frequency domain loop 
shaping requirements for the autopilot have been expressed as the performance 
functions with hard bound constraints. These performance functions are related with 
the disturbance rejection over the bandwidth of the feedback system, stability robustness 
to plant uncertainty and the shapings of the sensitivity and complementary sensitivity 
functions. These performance functions with constraints have been used to pose the 
frequency loop shaping problem as constrained optimization problem which has been 
solved using the Simulated Annealing (SA) technique. For getting the desired frequency 
loop shaping, the bounds on the performance functions are tuned during optimization 
process. The SA technique has been found suitable for achieving optimal solutions. 
The design methodology and algorithm developed here provide systematic synthesis 
procedure but the order of the autopilot turns out to be high. 
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In Chapter 4, the lateral (roll-yaw) autopilot has been synthesized. The autopilot 
has been designed for the system to track the stability axis roll rate while keeping 
the sideslip minimum. The design methodology and algorithm are based on the control 
law which provides input-output decoupling and thus facilitates independent designs 
for roll rate and yaw rate channels. The design methodology developed here, facilitates 
the stability axis roll rate command tracking with varying capabilities in yaw. The 
control law also facilitates the assignment of poles and zeros through the selection 
of polynomials. The polynomials have been selected for the desired time domain and 
frequency domain responses. The time domain and frequency domain specifications 
have been defined as hard bound constraints on the time domain and frequency 
domain performance functions. The time domain performance functions are related 
with the tracking error, overshoot, control input and sideslip whereas the frequency 
domain performance functions are related with the disturbance rejection, robust stability 
and, sensitivity and complementary sensitivity functions. The frequency domain 
performance functions have been defined at the output channels of the roll-yaw 
missile plant. These performance functions have been used to pose the design problem 
as the weighted sum constrained optimization problem which has been solved using 
the Simulated Annealing technique. The order of the autopilot is high in this design 
methodology also. 

In Chapter 5, the robustness analysis has been carried out using singular values 
for the designs in Chapters 2, 3 and 4. The uncertainty models considered are the 
neglected actuator dynamics, mismodeled actuator dynamics and the neglected sensor 
dynamics. These uncertainties have been modeled as multiplicative uncertainties. In 
the robustness analysis, the autopilot sensitivity due to the neglected actuator dynamics, 
mismodeled actuator dynamics and neglected sensor dynamics has been evaluated. 
The minimum actuator and sensor natural frequencies have been found out to guarantee 
the stability of the closed loop system. 

In summary, the contributions of this thesis are as follows : 

(1) A systematic design methodology and algorithm have been developed for 
the parametric synthesis of longitudinal autopilot for a BTT missile using direct 
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search methods of nonlinear optinnization. The autopilot has been designed taking 
one channel at a time. 

(2) A two degrees of freedom longitudinal autopilot has been designed in 
which normal acceleration command tracking and feedback properties have been 
independently achieved. A systematic design algorithm has been developed to achieve 
the desired command tracking using closed loop model matching concept, and to 
achieve the desired frequency domain specifications using constrained optimization. 
The constrained optimization problem for the frequency response shaping has been 
solved using Simulated Annealing technique. 

(3) A systematic design methodology and algorithm have been developed for 
the synthesis of input-output decoupled lateral (roll-yaw) autopilot for a BTT missile 
to track the stability axis roll rate command with minimum sideslip. The design 
methodology proposed, facilitates independent designs for the roll and yaw channels. 
In the design methodology proposed, various time domain and frequency domain 
specifications have been achieved by formulating the design problem as constrained 
optimization problem which has been solved using Simulated Annealing technique. 

(4) Singular value based analysis has been used to investigate the stability 
robustness against the neglected actuator dynamics, mismodeled actuator dynamics 
and the neglected sensor dynamics. The minimum actuator and sensor natural frequencies 
have been found out to guarantee stability of the closed loop system. 


6.3 Suggestions for Further Work 

In continuation of the research work carried out and reported in this thesis, the 
following problems need further investigations. 

(i) In all the three designs in Chapters 2, 3 and 4, the controller orders have 
come out to be high. The controller order reduction techniques may be applied to 
get the reduced order controllers. 

(ii) In Chapters 2 and 3, the loop shaping has been carried out at the input 
of the plant whereas in Chapter 4, the loop shaping has been earned out at the 
output of the plant. The design algorithms may be modified to incorporate the loop 
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shaping at both the input and output channels of the plant to improve the time 
domain and frequency domain responses. 

(iii) All the designs in the thesis have been carried out at one operating point. 
The feasibility of gain scheduling methods can be studied in the reported framework 
of the design. 
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Appendix 1 

Missile Dynamics 


The dynamics of a missile is described by a set of three force equations along 
x,y and z directions and three moment equations about x, y and z axis. This forms 
a set of nonlinear coupled equations. For a BTT missile, as shown in Fig. 1.1 of 
Chapter 1, the coupling is more prominent due to asymmetric nature of the airframe 
configuration. The nonlinear equations are linearized in order to form a linear time 
invariant (LTI) model. The linearized model so obtained is used for the design of 
LTI autopilot. 

Linearized missile dynamics for a kind of BTT missile shown in Fig. 1.1 is 
taken from [7] and is described in the following. 

Fig. 1.1 displays a highly asymmetric BTT missile configuration. Aerodynamic 
analysis of this airframe configuration shows that very strong roll-yaw coupling is 
present. Large roll rates are induced by sideslip angle created primarily by the 
asymmetry of the vehicle. For this reason, the linearized pitch dynamics are separated 
from linear coupled roll-yaw dynamics. The resulting autopilot design consists of a 
longitudinal pitch acceleration command autopilot and a lateral roll-yaw autopilot in 
which a roll rate about the velocity vector is commanded. 

The important variables defined in the Fig.1.1 are (p,q,r) which are body axis 
roll, pitch and yaw angular rates, respectively; (({>, 0, t]/) are integrals of (p,q,r); a is 
the angle of attack; and |3 is the sideslip angle. The LTI flight control design equations 
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arc foirncd by lincsrizing the nonlinear flerodynsmic eejuations of the motion about 
a trim condition. 

The pitch LTI flight control equations are : 


d = Za-oc + Z5-8^ + ^ 

(Al.l) 

q = M ■ a + • 5 

^ a 5 e 

(A1.2) 

5 = -2^0)5 -co^5 +0)^5 

€ e e ec 

(A1.3) 


The measurements that are available from the inertial measurement unit are 
normal acceleration Aj = VZa • cx + VZs • 8e and pitch rate g (rad/s). ■ The scalar 
control input u = dec (rad) is the elevator fin angle command. V, ^ and to are the 
missile velocity, actuator (elevator) damping and natural frequency respectively. 

The variables Za, Zg, Ma and Ms are the dimensional aerodynamic stability 
derivatives derived from aerodynamic measurements (wind tunnel) of lift and pitching 
moment. 

The LTI roll-yaw flight control equations are 


P = p sin (a) - r cos (a) + Kp • P + Kga ■ 5a + Y^r ■ 5r 

(A1.4) 

p = L. • p + L. -5 + L, -5 

^ p ^ oa a or r 

(A1.5) 

'r = N.-^ + N^ -5 +N^ -5 

P oa a or r 

(A1.6) 

5 = -2^co-5’ -co"5 +(0^5 

a ^ a a ac 

(AL7) 

5 = -2^co5'- 0)^5 + (0^5 

r ^ r r rc 

(A1.8) 


The variables Fp, YSa, YSr, Lp, LSa, LSr, iVp, NSa and NSr are dimensional 
aerodynamic stability derivatives derived from aerodynamic measurements of lateral 
side force and, roll and yaw moments. Equations (A1.7) and (A1.8) describe the 
actuator dynamics for aileron and mdder fins, da and 5/- are aileron and mdder fin 
deflections respectively. 

The stability axis coordinate system is defined as the transformation of body 
axes to the stability axes using angle of attack (cc). Wind axis coordinate system is 
defined as the transformation of stability axes to wind axes using sideslip p. The 
missile velocity vector is pointed along the x-axis of the wind coordinate system, as 
shown in Fig. 1.1. In a BTT autopilot, induced sideslip is minimized by rolling the 
missile about the velocity vector. Assuming that P is very small, this is a rotation 
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about the x-axis in the stability axis coordinate system. Thus the commanded roll 

rate is the stability axis roll rate ps, ie. 

Ps - p cos (a) + r sin (a) (A1 .9) 



Appendix 2 

Parametrization of Stabilizing 
Controllers with Two Degrees 
of Freedom 


A2(a) Servo Problem And Solvability Condition 


The system shown is Fig. 3.1 is said to be a servo system if the following 
two conditions are satisfied [42] : 


(i) Internal stability, meaning BIBO stability of the feedback system , ie. 


_C(/+GC)~‘ (Z + CG)’ _ 


(A2.1) 


We say that the pair H{G, C) is stable and that C stabilizes G if equation 
(A2.1) holds ie. all the transfer function matrices of matrix H{G, Q are stable rational 
matrices. 


(ii) Output regulation, meaning BIBO stability of the system with input ro and 
output e, ie. 

Hero is) = [ / + Mil + PC2)~ ’ PCi P * Z? 6 R_ (A2.2) 

The controller achieving these two properties is called the servo controller and 
the servo problem (P,M,R) is said to be solvable if such a controller exists. 


Let 


^ >r I P 
Pi } q-p 


ND"'- 


Ni 

N2 




(A2.3) 
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p q-p 


{A2A) 


be a right coprime factorization (r.c.O and a left coprime factorization (l.c.f.) 
respectively, of the plant and let /? = oc/ = \j/(t) ^ 7 be a coprime factorization with 


a , X 

(j)(5) 0 for s & C^ = |5:Re5<0} .The condition for solvability of the servo 

problem {P.M.R) given in [42] as theorem 1 is reproduced below as Lemma without 
proof : 


Lemma A2.1 

The servo problem {P,M,R) is solvable if and only if 

Ni and (^Ip are left coprime (A2.5) 

or equivalently. 

[N Dil and ([) Iq are left coprime (A2.6) 


A 2(b) The Class of Stabilizing Controllers 


This section gives two lemmas drawn from [42] for the servo system shown 
in Fig. 3.2 which contains the internal model I / ^{s) for the reference command at 
error channel. These two lemmas stated below without proof give the coprime 
factorization of Ge{s) and 


Pyis)/^s) 

P2{S) 


(A2.7) 


Lemma A2.2 

Suppose that l.c.f. of P{s) is represented by equation (A2.4) and that the 
condition (A2.5) or (A2.6) holds, then 

'<)/ 'o' 

Ge = be ^ Ne = 

Q b N (A2.8) 
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and 

F 4 , = iv<j, = [ 4 >Di D 2 j’ iv are l.c.f. 


(A2.9) 


Lemma A2.3 


Suppose that an r.c.f of P(s) is represented by (A 2.3) and that Ni(s)/i^{s) has 
an r.c.f. of 

Ai/(j) = i^D"^ (A2.10) 

then 

" - 


G = N D ^ = <t>^ (!>£))' 


and 




- - 

P.=N.D.'^= . (D b) ^ 

*•’ Nb 

2 


(A2.11) 


(A2.12) 


are r.c.f of and under the assumption of (A2.5) or (A2.6). 


The above two lemmas A2.2 and A2.3 give the possible classes of servo 
controllers and attainable characteristics stated in Lemma 3.1 of Chapter 3. The proofs 
of all above stated lemmas are available in [42]. 



Appendix 3 


Fundamental Relations of 
Sensitivity And Complementary 
Sensitivity Functions 

Fig.3.3, which is redrawn from Fig.3.1 in the standard form [32], including 
disturbance signal vector d(s) and measurement noise vector m{s). Other notations 
are defined in Section 3.2.2 with reference to Fig.3.1 of Chapter 3. 

From the Fig.3.3 we see that 

y{s) = d{s) + G{s) C{s) [r{s)-m{s) (A3.1) 

Hence, 

[/ + G(s) C(5)] y{s) = dis) + G{s) C{s) [/-(i) - ot(5)] (A3. 2) 

or 

3'(^) = [/ + G(s) C(^)]- ‘ {d{s) + G(s) C(s) [rfj) - (A3.3) 

Output sensitivity and complementary sensitivity functions are defined as 

So (s) = [/ + G{s) C(i)]" ‘ (A3.4) 

and 

To(s) = [/ + G{s) C(s)r ' G(j) C{s) , (A3.5) 

respectively. ( 

Then equation (A3.3) can be written as 

= Sois) d{s) + To(s) (A3.6) 

From the Fig.3.3 we see that 

«(^) = Cis) [ris) - mis) - d(s) - G{s) uis)] . 


(A3.7) 
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Hence, 

[/ + C{s) G(5)] uis) = C(s) [r(i) - m(s) - d{s)] . (A3 . 8) 

The input sensitivity function, Si{s), and input complementary sensitivity function 
Ti{s), are defined as 

Si{s) = [l + C{s)G{s)r^ (A3.9) 

and 

Ti{s) = C{s) G{s) [I + C(^) G{s)r ' (A3. 10) 

Then equation (A 3.8) can be written as 

u{s) = Si{s) C{s) [ris) - m{s) - d{s)] (A3. 11) 

From equations (A3 .4), (A3. 5), (A3. 9) and (A3. 10), the following two relations 
are obtained : 

Tiis) + Si{s) = I (A3. 12) 

Tq(s) + So(s) = ! (A3. 13) 

Dropping the suffix, in general, we can write 

T(s) + Sis) = I (A3. 14) 


From the equations (A3. 3) and (A3. 11) we see that for disturbance attenuation 
and for the output to be insensitive to the measurement noises, it is needed that 
both S{s) and T(s) should be small in some sense. However, equation (A3. 14) implies 
that both S{s) and T{s) cannot be made small simultaneously, if T(s) is made nearly 
zero, then S(s) becomes nearly identity. On the other hand, if S{s) is made nearly 
zero, then T(s) becomes almost identity. Thus we have an unavoidable trade off 
between the attenuating disturbances and filtering out the measurement noises. This 
trade off is one of the factors which make the feedbeck design difficult. The equation 
(A3. 14) is called fundamental relation between sensitivity and complementary 
sensitivity functions in feedback control. 



Appendix 4 

Simulated Annealing 


A4(a) Introduction 

The SA optimization algorithm can be considered analogous to the physical 
process by which a material changes state while minimising its energy [52]. A slow 
careful cooling brings the material to a highly ordered crystalline state of lowest 
energy. A rapid cooling instead yields defects and glass like intrusions inside the 
material. 

From optimization point of view, an iterative search accepting only new points 
with lower function values is like rapidly quenching a physical system at zero 
temperature. It is very likely to get stuck in a metastable local minimum. On the 
contrary, SA permits uphill moves under the control of a temperature parameter. At 
higher temperature only the gross behaviour of the cost function is relevant to the 
search. As temperature decreases finer details can be developed to get a good final 
point. 

The SA method is based on the random evaluations in such a way that transitions 
out of a local minimum are possible. This method is able to discriminate between 
"gross behaviour" of the function and finer "wrinkles". First, it reaches an area in 
the function domain where a global minimum should be present, following the gross 
behaviour irrespectively of small local minima found on the way. It then develops 
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finer details, finding a good, near optimal local minima, if not the global minimum 
itself. The details of SA algorithm in the following sections are drawn from [52], 

A4(b) Method 

Let X be a vector in i?" and (,vi , .X2 , .... , x„) its components. Let f(X) be the 
function to minimize and let a\ < x\ < b\ , ...., < Xn < bn be its n variables each 

ranging in a finite continuous interval. The function f does not need to be bounded 

The SA algorithm is schematically shown in Fig. A4.1. It proceeds iteratively 
: Starting from a given point Xq it generates a succession of points : Xq, X\, ...., Xi, 
... tending to the global minimum of the cost function. New candidate points are 
generated around the current point Xi applying random moves along each coordinate 
direction, in turn. The new coordinate values are uniformly distributed in intervals 
centered round the corresponding coordinate of Xi. Half the size of these intervals 
along each coordinate is recorded in the step vector V. If the point falls outside the 
definition domain of f, a new point is randomly generated until a point belonging 
to the definition domain is found. A candidate point X is accepted or rejected 
according to Metropolis Criterion [58] : 

If A/ < 0, then accept the new point : Xi+i = X' else accept the new point 

with probability: 

p (A/) = exp (- Af/T) , 

where A/ = f{X ') -f {Xi) and T is a parameter called temperature. 

At a fixed value of T the succession of points Xq , Xi ,...., X/ , .... is not 
downhill, except when T=0. For values of T large compared to the mean value of 
|/(X/,)-/(X^) I (X/i and Xk are points randomly chosen inside the definition domain 
of f) almost all new points are accepted and the succession is a random sampling 
of f. 

The SA algorithm starts at some "high" temperature To given by the user. A 
sequence of points is then generated until a sort of "equilibrium" is approached ; 
that is a sequence of points Xi whose average value of f reaches a stable value as 
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Initialize parameters 

4 

Perform a cycle of random moves, each along 
one co-ordinate direction. Accept or reject 
each point according to the Metropolis criterion. 
Record the optimum point reached so far. 



Adjust step vector V. 
Reset No. of Cycles to 0. 



Reduce temperature. 

Reset No. of adjustments to 0. 

Set current point to the optimum. 



( END ) 


A /f t - r» A 
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i increases. During this phase the step vector Vm is periodically adjusted to better 
follow the function behaviour. The best point reached is recorded as Hopt. 

After thermal equilibration, the temperature T is reduced and a new sequence 
of moves is made starting from 'Kopt, until the thermal equilibrium is reached again, 
and so on. The process is stopped at a temperature low enough that no more useful 
improvement can be expected, according to a stopping criterion that will be described 
later. 

The detailed description of the algorithm is given below : 

Step 0 (Initialization) 

Choose 

A starting point Xo . 

A starting step vector Vq. 

A starting temperature Tq. 

A termination criterion 8 and a number of successive temperature reductions to 
test for termination Ne .A test for step variation Ns and a varying criterion c. 

A test for temperature reduction Nt and a reduction coefficient rj. 

Set i, j, m and k to 0, i is the index denoting successive points, j denotes 
successive cycles along every direction, m describes successive step adjustments, 
and k covers successive temperature reductions. 

Set /z to 1. h is the index denoting the direction along which the trial point 
is generated, starting from the last accepted point. 

Compute /o = /(Xo) 

Set Xopt = Xo , fopt = fo ■ 

Set TZu = 0 , M = 1, ... , /I . 

Stifu=fo,u = . 

Step 1 

Starting from the point X/, generate a random point X' along the direction h 
as 

X = Xi + rVmh Sh 


151 


where r is a random number generated in the range [-1,1] by a pseudo random 
number generator; eh is the vector of the hth co-ordinate direction; and V,nh is the 
component of the step vector Vm along the same direction. 

Step 2 

If the hth coordinate of X' lies outside the definition domain of f, that is, if 
J// < ah or Xh' > bh , then return to Step 1. 

Step 3 

Compute f' - f{X') . 

If f' < fi, then accept the new point ; 

set Xi + 1 = X ', 

set fi+l = f\ 

add 1 to f , 

add 1 to nh ; 

If /' < fopt, then set 
Xopt = X' 

fopt = f' ■ 

endif; 

else if' >fi) accept or reject the point with acceptance probability p (Metropolis 
move). 



In practice, a pseudorandom number p' is generated in the range [0,1] and is 
compared with p. If p' < p, the point is accepted, otherwise it is rejected. 

In the case of acceptance 
set Xf+i = X' , 
set/j+i = /' , 
add 1 to I , 

add 1 to Till. . 

Step 4 

Add 1 to h. 


I 
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If h < n, then go to Step 1 ; 
else set h to 1 and add 1 to 7 . 

Step 5 

If j < Ns then go to step 1 ; 
else update the step vector V,n\ 

for each direction u the new step vector component Yu ' is 
^ ?iu/Ns - 0.6 ^ 

Vu' = Vnm 1 + C„ Q ^ if > 0.6 Ns , 

T/ / -r r\ A XT 

Vu — ^ . if riu < 0.4 A^5 , 

, ^ ^ 0.4 - ?iu/Ns 

^ 0.4 

vm ' = Vmu otherwise 

Set V,„+i = V', 
set j to 0, 

set flu to 0, M = 1 , .... , n, 
add 1 to m . 

The aim of these variations in step length is to maintain the average percentage 
of accepted moves at about one-half of the total number of moves. The Cu parameter 
controls the step variaton along each uth direction. 

Step 6 

If m < Nt, then go to step 1; 

else, it is time to reduce the temperature Tk : 

set Tk+ 1 = nr - Tk, 

set fk =fi , 
add 1 to ^ , 
set m to 0. 

It is worth noting that a temperature reduction occurs every Ns • Nj cycles of 
moves along every direction and after Nj step adjustments. 
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Step 7 (termination criterion) 

If : 

fk-fk-u ^ e, « = 1, 
and 

fii ~fopi ^ e 

then stop the search; 
else : 

add 1 to i. 

Set Xi = Xopt 1 
set fi = fopt ■ 

Go to Step 1. 



